MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
adjoint_advection_diffusion.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // ----------------------------------------------------------
13 // Advection-Diffusion Miniapp: Parallel MFEM CVODES Example
14 // ----------------------------------------------------------
15 //
16 // Compile with: make adjoint_advection_diffusion
17 //
18 // Sample runs: adjoint_advection_diffusion -dt 0.01 -tf 2.5
19 // adjoint_advection_diffusion -dt 0.005
20 //
21 // Description: This example is a port of cvodes/parallel/cvsAdvDiff_ASAp_non_p
22 // example that is part of SUNDIALS. The goal is to demonstrate
23 // how to use the adjoint SUNDIALS CVODES interface with MFEM.
24 // Below is an excerpt description from the aforementioned file.
25 //
26 // Example problem:
27 //
28 // The following is a simple example problem, with the program for its solution
29 // by CVODE. The problem is the semi-discrete form of the advection-diffusion
30 // equation in 1-D:
31 //
32 // du/dt = p1 * d^2u / dx^2 + p2 * du / dx
33 //
34 // on the interval 0 <= x <= 2, and the time interval 0 <= t <= 5. Homogeneous
35 // Dirichlet boundary conditions are posed, and the initial condition is:
36 //
37 // u(x,t=0) = x(2-x)exp(2x).
38 //
39 // The nominal values of the two parameters are: p1=1.0, p2=0.5. The PDE is
40 // discretized on a uniform grid of size MX+2 with central differencing, and
41 // with boundary values eliminated, leaving an ODE system of size NEQ = MX.
42 //
43 // The program solves the problem with the option for nonstiff systems: ADAMS
44 // method and functional iteration. It uses scalar relative and absolute
45 // tolerances. In addition to the solution, sensitivities with respect to p1 and
46 // p2 as well as with respect to initial conditions are computed for the
47 // quantity:
48 //
49 // g(t, u, p) = int_x u(x,t) at t = 5
50 //
51 // These sensitivities are obtained by solving the adjoint system:
52 //
53 // dv/dt = -p1 * d^2 v / dx^2 + p2 * dv / dx
54 //
55 // with homogeneous Dirichlet boundary conditions and the final condition:
56 //
57 // v(x,t=5) = 1.0
58 //
59 // Then, v(x, t=0) represents the sensitivity of g(5) with respect to u(x, t=0)
60 // and the gradient of g(5) with respect to p1, p2 is
61 //
62 // (dg/dp)^T = [ int_t int_x (v * d^2u / dx^2) dx dt ]
63 // [ int_t int_x (v * du / dx) dx dt ]
64 //
65 // This version uses MPI for user routines.
66 // Execute with number of processors = N, with 1 <= N <= MX.
67 
68 #include "mfem.hpp"
69 #include <fstream>
70 #include <iostream>
71 #include <algorithm>
72 
73 #ifndef MFEM_USE_SUNDIALS
74 #error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
75 #endif
76 
77 using namespace std;
78 using namespace mfem;
79 
80 // Implement the adjoint rate equations in AdvDiffSUNDIALS
81 class AdvDiffSUNDIALS : public TimeDependentAdjointOperator
82 {
83 public:
84  AdvDiffSUNDIALS(int ydot_dim, int ybdot_dim, Vector p,
85  ParFiniteElementSpace *fes, Array<int> & ess_tdof) :
86  TimeDependentAdjointOperator(ydot_dim, ybdot_dim),
87  p_(p),
88  ess_tdof_list(ess_tdof),
89  pfes(fes),
90  Mf(NULL),
91  M_solver(fes->GetComm())
92  {
93  int skip_zeros = 0;
94 
95  cout << "Essential tdofs: " << endl;
96  ess_tdof_list.Print();
97 
98  m = new ParBilinearForm(pfes);
99  m->AddDomainIntegrator(new MassIntegrator());
100  m->Assemble(skip_zeros);
101  m->Finalize(skip_zeros);
102 
103  // Define coefficients
104  mp0 = new ConstantCoefficient(-p_[0]);
105  p0 = new ConstantCoefficient(p_[0]);
106  Vector p2vec(fes->GetParMesh()->SpaceDimension());
107  p2vec = p_[1];
108  p2 = new VectorConstantCoefficient(p2vec);
109 
110  k = new ParBilinearForm(pfes);
111  k->AddDomainIntegrator(new DiffusionIntegrator(*mp0));
112  k->AddDomainIntegrator(new ConvectionIntegrator(*p2));
113  k->Assemble(skip_zeros);
114  k->Finalize(skip_zeros);
115 
116  k1 = new ParBilinearForm(pfes);
117  k1->AddDomainIntegrator(new DiffusionIntegrator(*p0));
118  k1->AddDomainIntegrator(new ConvectionIntegrator(*p2));
119  k1->Assemble(skip_zeros);
120  k1->Finalize(skip_zeros);
121 
122  M = m->ParallelAssemble();
123  HypreParMatrix *temp = M->EliminateRowsCols(ess_tdof_list);
124  delete temp;
125 
126  K = k->ParallelAssemble();
127  temp = K->EliminateRowsCols(ess_tdof_list);
128  delete temp;
129 
130  K_adj = k1->ParallelAssemble();
131  temp = K_adj->EliminateRowsCols(ess_tdof_list);
132  delete temp;
133 
134  M_prec.SetType(HypreSmoother::Jacobi);
135  M_solver.SetPreconditioner(M_prec);
136  M_solver.SetOperator(*M);
137 
138  M_solver.SetRelTol(1e-14);
139  M_solver.SetAbsTol(0.0);
140  M_solver.SetMaxIter(1000);
141  M_solver.SetPrintLevel(0);
142 
143  }
144 
145  virtual void Mult(const Vector &x, Vector &y) const;
146 
147  virtual void AdjointRateMult(const Vector &y, Vector &yB,
148  Vector &yBdot) const;
149 
150  virtual int SUNImplicitSetup(const Vector &y,
151  const Vector &fy, int jok, int *jcur,
152  double gamma);
153 
154  virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
155 
156  virtual void QuadratureSensitivityMult(const Vector &y, const Vector &yB,
157  Vector &qbdot) const;
158 
159  ~AdvDiffSUNDIALS()
160  {
161  delete m;
162  delete k;
163  delete k1;
164  delete M;
165  delete K;
166  delete K_adj;
167  delete Mf;
168  delete p0;
169  delete mp0;
170  delete p2;
171  }
172 
173 protected:
174  Vector p_;
175  Array<int> ess_tdof_list;
176  ParFiniteElementSpace *pfes;
177 
178  // Internal matrices
179  ParBilinearForm *m;
180  ParBilinearForm *k;
181  ParBilinearForm *k1;
182 
183  HypreParMatrix *M;
184  HypreParMatrix *K;
185  HypreParMatrix *K_adj;
186 
187  HypreParMatrix *Mf;
188 
189  CGSolver M_solver;
190  HypreSmoother M_prec;
192  ConstantCoefficient *mp0;
194 };
195 
196 // Initial conditions for the problem
197 double u_init(const Vector &x)
198 {
199  return x[0]*(2. - x[0])*exp(2.*x[0]);
200 }
201 
202 int main(int argc, char *argv[])
203 {
204  // Initialize MPI and HYPRE.
205  Mpi::Init(argc, argv);
206  int myid = Mpi::WorldRank();
207  Hypre::Init();
208 
209  // Parse command-line options.
210  int ser_ref_levels = 0;
211  int par_ref_levels = 0;
212  double t_final = 2.5;
213  double dt = 0.01;
214  int mx = 20;
215  bool step_mode = true;
216 
217  int precision = 8;
218  cout.precision(precision);
219 
220  OptionsParser args(argc, argv);
221 
222  args.AddOption(&mx, "-m", "--mx", "The number of mesh elements in the x-dir");
223  args.AddOption(&ser_ref_levels, "-r", "--refine",
224  "Number of times to refine the mesh uniformly.");
225  args.AddOption(&step_mode, "-a", "--adams", "-no-a","--no-adams",
226  "A switch to toggle between CV_ADAMS, and CV_BDF stepping modes");
227  args.AddOption(&t_final, "-tf", "--t-final",
228  "Final time; start time is 0.");
229  args.AddOption(&dt, "-dt", "--time-step",
230  "Time step.");
231 
232  args.Parse();
233  if (!args.Good())
234  {
235  if (myid == 0)
236  {
237  args.PrintUsage(cout);
238  }
239  return 1;
240  }
241  if (myid == 0)
242  {
243  args.PrintOptions(cout);
244  }
245 
246  // Create a small 1D mesh with a length of 2. This mesh corresponds with the
247  // cvsAdvDiff_ASA_p_non_p example.
248  Mesh mesh = Mesh::MakeCartesian1D(mx+1, 2.);
249 
250  // Refine the mesh to increase the resolution. In this example we do
251  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
252  // command-line parameter. If the mesh is of NURBS type, we convert it to
253  // a (piecewise-polynomial) high-order mesh.
254  for (int lev = 0; lev < ser_ref_levels; lev++)
255  {
256  mesh.UniformRefinement();
257  }
258 
259  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
260  mesh.Clear();
261  for (int lev = 0; lev < par_ref_levels; lev++)
262  {
263  pmesh->UniformRefinement();
264  }
265 
266  // Finite Element Spaces
267  H1_FECollection fec(1, pmesh->SpaceDimension());
268  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
269 
270  HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
271  if (myid == 0)
272  {
273  cout << "Number of unknowns: " << global_vSize << endl;
274  }
275 
276  // Set up material properties, and primal and adjoint variables
277  // p are the fixed material properties
278  Vector p(2);
279  p[0] = 1.0;
280  p[1] = 0.5;
281 
282  // U is the size of the solution/primal vector
283  // Set U with the initial conditions
284  ParGridFunction u(fes);
286  u.ProjectCoefficient(u0);
287 
288  cout << "Init u: " << endl;
289  u.Print();
290 
291  // TimeDependentOperators need to be TrueDOF Size
292  HypreParVector *U = u.GetTrueDofs();
293 
294  // Get boundary conditions
295  Array<int> ess_tdof_list;
296  Array<int> essential_attr(pmesh->bdr_attributes.Size());
297  essential_attr[0] = 1;
298  essential_attr[1] = 1;
299  fes->GetEssentialTrueDofs(essential_attr, ess_tdof_list);
300 
301  // Setup the TimeDependentAdjointOperator and the CVODESSolver
302  AdvDiffSUNDIALS adv(U->Size(), U->Size(), p, fes, ess_tdof_list);
303 
304  // Set the initial time to the TimeDependentAdjointOperator
305  double t = 0.0;
306  adv.SetTime(t);
307 
308  // Create the CVODES solver corresponding to the selected step method
309  CVODESSolver *cvodes = new CVODESSolver(fes->GetComm(),
310  step_mode ? CV_ADAMS : CV_BDF);
311  cvodes->Init(adv);
312  cvodes->UseSundialsLinearSolver();
313  cvodes->SetMaxNSteps(5000);
314 
315  // Relative and absolute tolerances for CVODES
316  double reltol = 1e-8, abstol = 1e-6;
317  cvodes->SetSStolerances(reltol, abstol);
318 
319  // Initialize adjoint problem settings
320  int checkpoint_steps = 50; // steps between checkpoints
321  cvodes->InitAdjointSolve(checkpoint_steps, CV_HERMITE);
322 
323  // Perform time-integration for the problem (looping over the time
324  // iterations, ti, with a time-step dt).
325  bool done = false;
326  for (int ti = 0; !done; )
327  {
328  double dt_real = max(dt, t_final - t);
329  cvodes->Step(*U, t, dt_real);
330  ti++;
331 
332  done = (t >= t_final - 1e-8*dt);
333 
334  if (done && myid == 0)
335  {
336  cvodes->PrintInfo();
337  }
338  }
339 
340  u = *U;
341  if (myid == 0)
342  {
343  cout << "Final Solution: " << t << endl;
344  }
345 
346  cout << "u (" << myid << "):" << endl;
347  u.Print();
348  cout << flush;
349  MPI_Barrier(MPI_COMM_WORLD);
350 
351  // Calculate the quadrature int_x u dx at t = 5
352  // Since it's only a spatial quadrature we evaluate it at t=5
353  ParLinearForm obj(fes);
354  ConstantCoefficient one(1.0);
356  obj.Assemble();
357 
358  double g = obj(u);
359  if (myid == 0)
360  {
361  cout << "g: " << g << endl;
362  }
363 
364  // Solve the adjoint problem. v is the adjoint solution
365  ParGridFunction v(fes);
366  v = 1.;
367  v.SetSubVector(ess_tdof_list, 0.0);
368  HypreParVector *V = v.GetTrueDofs();
369 
370  // Initialize quadrature sensitivity values to zero
371  Vector qBdot(p.Size());
372  qBdot = 0.;
373 
374  t = t_final;
375  cvodes->InitB(adv);
376  cvodes->InitQuadIntegrationB(qBdot, 1.e-6, 1.e-6);
377  cvodes->UseSundialsLinearSolverB();
378  cvodes->SetSStolerancesB(reltol, abstol);
379 
380  // Results at time TBout1
381  double dt_real = max(dt, t);
382  cvodes->StepB(*V, t, dt_real);
383  if (myid == 0)
384  {
385  cout << "t: " << t << endl;
386  }
387 
388  cout << "v (" << myid << "):" << endl;
389  V->Vector::Print();
390  cout << flush;
391  MPI_Barrier(MPI_COMM_WORLD);
392 
393  // Evaluate the sensitivity
394  cvodes->EvalQuadIntegrationB(t, qBdot);
395 
396  MPI_Barrier(MPI_COMM_WORLD);
397  if (myid == 0)
398  {
399  cout << "sensitivity:" << endl;
400  qBdot.Print();
401  }
402 
403  // Free the used memory.
404  delete fes;
405  delete pmesh;
406  delete U;
407  delete V;
408  delete cvodes;
409 
410  return 0;
411 }
412 
413 // AdvDiff rate equation
414 void AdvDiffSUNDIALS::Mult(const Vector &x, Vector &y) const
415 {
416  Vector z(x.Size());
417  Vector x1(x);
418 
419  // Set boundary conditions to zero
420  x1.SetSubVector(ess_tdof_list, 0.0);
421 
422  K->Mult(x1, z);
423 
424  y = 0.;
425  M_solver.Mult(z, y);
426 }
427 
428 // AdvDiff Rate equation setup
429 int AdvDiffSUNDIALS::SUNImplicitSetup(const Vector &y,
430  const Vector &fy,
431  int jok, int *jcur, double gamma)
432 {
433  // Mf = M(I - gamma J) = M - gamma * M * J
434  // J = df/dy => K
435  *jcur = 1; // We've updated the jacobian
436 
437  delete Mf;
438  Mf = Add(1., *M, -gamma, *K);
439  HypreParMatrix *temp = Mf->EliminateRowsCols(ess_tdof_list);
440  delete temp;
441  return 0;
442 }
443 
444 // AdvDiff Rate equation solve
445 int AdvDiffSUNDIALS::SUNImplicitSolve(const Vector &b, Vector &x, double tol)
446 {
447  Vector z(b.Size());
448  M->Mult(b,z);
449 
450  CGSolver solver(pfes->GetComm());
451  HypreSmoother prec;
452  prec.SetType(HypreSmoother::Jacobi);
453  solver.SetPreconditioner(prec);
454  solver.SetOperator(*Mf);
455  solver.SetRelTol(1E-14);
456  solver.SetMaxIter(1000);
457 
458  solver.Mult(z, x);
459 
460  return (0);
461 }
462 
463 // AdvDiff adjoint rate equation
464 void AdvDiffSUNDIALS::AdjointRateMult(const Vector &y, Vector & yB,
465  Vector &yBdot) const
466 {
467  Vector z(yB.Size());
468 
469  // Set boundary conditions to zero
470  yB.SetSubVector(ess_tdof_list, 0.0);
471  K_adj->Mult(yB, z);
472  M_solver.Mult(z, yBdot);
473 }
474 
475 // AdvDiff quadrature sensitivity rate equation
476 void AdvDiffSUNDIALS::QuadratureSensitivityMult(const Vector &y,
477  const Vector &yB,
478  Vector &qBdot) const
479 {
480  // Now we have both the adjoint, yB, and y, at the same point in time
481  // We calculate
482  /*
483  * to u(x, t=0) and the gradient of g(5) with respect to p1, p2 is
484  * (dg/dp)^T = [ int_t int_x (v * d^2u / dx^2) dx dt ]
485  * [ int_t int_x (v * du / dx) dx dt ]
486  */
487 
488  ParBilinearForm dp1(pfes);
489  ConstantCoefficient mone(-1.);
490  dp1.AddDomainIntegrator(new DiffusionIntegrator(mone));
491  dp1.Assemble();
492  dp1.Finalize();
493 
494  HypreParMatrix * dP1 = dp1.ParallelAssemble();
495  HypreParMatrix *temp = dP1->EliminateRowsCols(ess_tdof_list);
496  delete temp;
497 
498  Vector b1(y.Size());
499  dP1->Mult(y, b1);
500  delete dP1;
501 
502  ParBilinearForm dp2(pfes);
503  Vector p2vec(pfes->GetParMesh()->SpaceDimension()); p2vec = 1.;
504  VectorConstantCoefficient dp2_coef(p2vec);
505  dp2.AddDomainIntegrator(new ConvectionIntegrator(dp2_coef));
506  dp2.Assemble();
507  dp2.Finalize();
508 
509  HypreParMatrix * dP2 = dp2.ParallelAssemble();
510  temp = dP2->EliminateRowsCols(ess_tdof_list);
511  delete temp;
512 
513  Vector b2(y.Size());
514  dP2->Mult(y, b2);
515  delete dP2;
516 
517  double dp1_result = InnerProduct(pfes->GetComm(), yB, b1);
518  double dp2_result = InnerProduct(pfes->GetComm(), yB, b2);
519 
520  qBdot[0] = -dp1_result;
521  qBdot[1] = -dp2_result;
522 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:2240
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:560
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
Conjugate gradient method.
Definition: solvers.hpp:465
virtual void StepB(Vector &w, double &t, double &dt)
Solve one adjoint time step.
Definition: sundials.cpp:1124
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:698
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
Vector coefficient that is constant in space and time.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
FDualNumber< tbase > exp(const FDualNumber< tbase > &f)
exp([dual number])
Definition: fdual.hpp:499
virtual void Step(Vector &x, double &t, double &dt)
Definition: sundials.cpp:1096
void UseSundialsLinearSolverB()
Use built in SUNDIALS Newton solver.
Definition: sundials.cpp:969
void SetSStolerancesB(double reltol, double abstol)
Tolerance specification functions for the adjoint problem.
Definition: sundials.cpp:1016
void InitAdjointSolve(int steps, int interpolation)
Initialize Adjoint.
Definition: sundials.cpp:895
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1916
Class for parallel linear form.
Definition: plinearform.hpp:26
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB=1e-3, double abstolQB=1e-8)
Initialize Quadrature Integration (Adjoint)
Definition: sundials.cpp:922
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
Parallel smoothers in hypre.
Definition: hypre.hpp:896
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
int SpaceDimension() const
Definition: mesh.hpp:1000
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:148
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:711
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
HYPRE_Int HYPRE_BigInt
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:426
void InitB(TimeDependentAdjointOperator &f_)
Initialize the adjoint problem.
Definition: sundials.cpp:862
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
Class for parallel bilinear form.
RefCoord t[3]
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:904
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1703
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:742
void Init(TimeDependentAdjointOperator &f_)
Definition: sundials.cpp:857
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
double u_init(const Vector &x)
Class for parallel meshes.
Definition: pmesh.hpp:32
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:3224
int main()
void EvalQuadIntegrationB(double t, Vector &dG_dp)
Evaluate Quadrature solution.
Definition: sundials.cpp:838
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition: sundials.cpp:678
void SetMaxNSteps(int steps)
Set the maximum number of time steps.
Definition: sundials.cpp:722
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
alpha (q . grad u, v)