MFEM  v4.3.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-2021, 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.
205  int num_procs, myid;
206  MPI_Init(&argc, &argv);
207  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
208  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
209 
210  // Parse command-line options.
211  int ser_ref_levels = 0;
212  int par_ref_levels = 0;
213  double t_final = 2.5;
214  double dt = 0.01;
215  int mx = 20;
216  bool step_mode = true;
217 
218  int precision = 8;
219  cout.precision(precision);
220 
221  OptionsParser args(argc, argv);
222 
223  args.AddOption(&mx, "-m", "--mx", "The number of mesh elements in the x-dir");
224  args.AddOption(&ser_ref_levels, "-r", "--refine",
225  "Number of times to refine the mesh uniformly.");
226  args.AddOption(&step_mode, "-a", "--adams", "-no-a","--no-adams",
227  "A switch to toggle between CV_ADAMS, and CV_BDF stepping modes");
228  args.AddOption(&t_final, "-tf", "--t-final",
229  "Final time; start time is 0.");
230  args.AddOption(&dt, "-dt", "--time-step",
231  "Time step.");
232 
233  args.Parse();
234  if (!args.Good())
235  {
236  if (myid == 0)
237  {
238  args.PrintUsage(cout);
239  }
240  MPI_Finalize();
241  return 1;
242  }
243  if (myid == 0)
244  {
245  args.PrintOptions(cout);
246  }
247 
248  // Create a small 1D mesh with a length of 2. This mesh corresponds with the
249  // cvsAdvDiff_ASA_p_non_p example.
250  Mesh mesh = Mesh::MakeCartesian1D(mx+1, 2.);
251 
252  // Refine the mesh to increase the resolution. In this example we do
253  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
254  // command-line parameter. If the mesh is of NURBS type, we convert it to
255  // a (piecewise-polynomial) high-order mesh.
256  for (int lev = 0; lev < ser_ref_levels; lev++)
257  {
258  mesh.UniformRefinement();
259  }
260 
261  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
262  mesh.Clear();
263  for (int lev = 0; lev < par_ref_levels; lev++)
264  {
265  pmesh->UniformRefinement();
266  }
267 
268  // Finite Element Spaces
269  H1_FECollection fec(1, pmesh->SpaceDimension());
270  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
271 
272  HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
273  if (myid == 0)
274  {
275  cout << "Number of unknowns: " << global_vSize << endl;
276  }
277 
278  // Set up material properties, and primal and adjoint variables
279  // p are the fixed material properties
280  Vector p(2);
281  p[0] = 1.0;
282  p[1] = 0.5;
283 
284  // U is the size of the solution/primal vector
285  // Set U with the initial conditions
286  ParGridFunction u(fes);
288  u.ProjectCoefficient(u0);
289 
290  cout << "Init u: " << endl;
291  u.Print();
292 
293  // TimeDependentOperators need to be TrueDOF Size
294  HypreParVector *U = u.GetTrueDofs();
295 
296  // Get boundary conditions
297  Array<int> ess_tdof_list;
298  Array<int> essential_attr(pmesh->bdr_attributes.Size());
299  essential_attr[0] = 1;
300  essential_attr[1] = 1;
301  fes->GetEssentialTrueDofs(essential_attr, ess_tdof_list);
302 
303  // Setup the TimeDependentAdjointOperator and the CVODESSolver
304  AdvDiffSUNDIALS adv(U->Size(), U->Size(), p, fes, ess_tdof_list);
305 
306  // Set the initial time to the TimeDependentAdjointOperator
307  double t = 0.0;
308  adv.SetTime(t);
309 
310  // Create the CVODES solver corresponding to the selected step method
311  CVODESSolver *cvodes = new CVODESSolver(fes->GetComm(),
312  step_mode ? CV_ADAMS : CV_BDF);
313  cvodes->Init(adv);
314  cvodes->UseSundialsLinearSolver();
315  cvodes->SetMaxNSteps(5000);
316 
317  // Relative and absolute tolerances for CVODES
318  double reltol = 1e-8, abstol = 1e-6;
319  cvodes->SetSStolerances(reltol, abstol);
320 
321  // Initialize adjoint problem settings
322  int checkpoint_steps = 50; // steps between checkpoints
323  cvodes->InitAdjointSolve(checkpoint_steps, CV_HERMITE);
324 
325  // Perform time-integration for the problem (looping over the time
326  // iterations, ti, with a time-step dt).
327  bool done = false;
328  for (int ti = 0; !done; )
329  {
330  double dt_real = max(dt, t_final - t);
331  cvodes->Step(*U, t, dt_real);
332  ti++;
333 
334  done = (t >= t_final - 1e-8*dt);
335 
336  if (done && myid == 0)
337  {
338  cvodes->PrintInfo();
339  }
340  }
341 
342  u = *U;
343  if (myid == 0)
344  {
345  cout << "Final Solution: " << t << endl;
346  }
347 
348  cout << "u (" << myid << "):" << endl;
349  u.Print();
350  cout << flush;
351  MPI_Barrier(MPI_COMM_WORLD);
352 
353  // Calculate the quadrature int_x u dx at t = 5
354  // Since it's only a spatial quadrature we evaluate it at t=5
355  ParLinearForm obj(fes);
356  ConstantCoefficient one(1.0);
358  obj.Assemble();
359 
360  double g = obj(u);
361  if (myid == 0)
362  {
363  cout << "g: " << g << endl;
364  }
365 
366  // Solve the adjoint problem. v is the adjoint solution
367  ParGridFunction v(fes);
368  v = 1.;
369  v.SetSubVector(ess_tdof_list, 0.0);
370  HypreParVector *V = v.GetTrueDofs();
371 
372  // Initialize quadrature sensitivity values to zero
373  Vector qBdot(p.Size());
374  qBdot = 0.;
375 
376  t = t_final;
377  cvodes->InitB(adv);
378  cvodes->InitQuadIntegrationB(qBdot, 1.e-6, 1.e-6);
379  cvodes->UseSundialsLinearSolverB();
380  cvodes->SetSStolerancesB(reltol, abstol);
381 
382  // Results at time TBout1
383  double dt_real = max(dt, t);
384  cvodes->StepB(*V, t, dt_real);
385  if (myid == 0)
386  {
387  cout << "t: " << t << endl;
388  }
389 
390  cout << "v (" << myid << "):" << endl;
391  V->Vector::Print();
392  cout << flush;
393  MPI_Barrier(MPI_COMM_WORLD);
394 
395  // Evaluate the sensitivity
396  cvodes->EvalQuadIntegrationB(t, qBdot);
397 
398  MPI_Barrier(MPI_COMM_WORLD);
399  if (myid == 0)
400  {
401  cout << "sensitivity:" << endl;
402  qBdot.Print();
403  }
404 
405  // Free the used memory.
406  delete fes;
407  delete pmesh;
408  delete U;
409  delete V;
410  delete cvodes;
411 
412  MPI_Finalize();
413  return 0;
414 }
415 
416 // AdvDiff rate equation
417 void AdvDiffSUNDIALS::Mult(const Vector &x, Vector &y) const
418 {
419  Vector z(x.Size());
420  Vector x1(x);
421 
422  // Set boundary conditions to zero
423  x1.SetSubVector(ess_tdof_list, 0.0);
424 
425  K->Mult(x1, z);
426 
427  y = 0.;
428  M_solver.Mult(z, y);
429 }
430 
431 // AdvDiff Rate equation setup
432 int AdvDiffSUNDIALS::SUNImplicitSetup(const Vector &y,
433  const Vector &fy,
434  int jok, int *jcur, double gamma)
435 {
436  // Mf = M(I - gamma J) = M - gamma * M * J
437  // J = df/dy => K
438  *jcur = 1; // We've updated the jacobian
439 
440  delete Mf;
441  Mf = Add(1., *M, -gamma, *K);
442  HypreParMatrix *temp = Mf->EliminateRowsCols(ess_tdof_list);
443  delete temp;
444  return 0;
445 }
446 
447 // AdvDiff Rate equation solve
448 int AdvDiffSUNDIALS::SUNImplicitSolve(const Vector &b, Vector &x, double tol)
449 {
450  Vector z(b.Size());
451  M->Mult(b,z);
452 
453  CGSolver solver(pfes->GetComm());
454  HypreSmoother prec;
455  prec.SetType(HypreSmoother::Jacobi);
456  solver.SetPreconditioner(prec);
457  solver.SetOperator(*Mf);
458  solver.SetRelTol(1E-14);
459  solver.SetMaxIter(1000);
460 
461  solver.Mult(z, x);
462 
463  return (0);
464 }
465 
466 // AdvDiff adjoint rate equation
467 void AdvDiffSUNDIALS::AdjointRateMult(const Vector &y, Vector & yB,
468  Vector &yBdot) const
469 {
470  Vector z(yB.Size());
471 
472  // Set boundary conditions to zero
473  yB.SetSubVector(ess_tdof_list, 0.0);
474  K_adj->Mult(yB, z);
475  M_solver.Mult(z, yBdot);
476 }
477 
478 // AdvDiff quadrature sensitivity rate equation
479 void AdvDiffSUNDIALS::QuadratureSensitivityMult(const Vector &y,
480  const Vector &yB,
481  Vector &qBdot) const
482 {
483  // Now we have both the adjoint, yB, and y, at the same point in time
484  // We calculate
485  /*
486  * to u(x, t=0) and the gradient of g(5) with respect to p1, p2 is
487  * (dg/dp)^T = [ int_t int_x (v * d^2u / dx^2) dx dt ]
488  * [ int_t int_x (v * du / dx) dx dt ]
489  */
490 
491  ParBilinearForm dp1(pfes);
492  ConstantCoefficient mone(-1.);
493  dp1.AddDomainIntegrator(new DiffusionIntegrator(mone));
494  dp1.Assemble();
495  dp1.Finalize();
496 
497  HypreParMatrix * dP1 = dp1.ParallelAssemble();
498  HypreParMatrix *temp = dP1->EliminateRowsCols(ess_tdof_list);
499  delete temp;
500 
501  Vector b1(y.Size());
502  dP1->Mult(y, b1);
503  delete dP1;
504 
505  ParBilinearForm dp2(pfes);
506  Vector p2vec(pfes->GetParMesh()->SpaceDimension()); p2vec = 1.;
507  VectorConstantCoefficient dp2_coef(p2vec);
508  dp2.AddDomainIntegrator(new ConvectionIntegrator(dp2_coef));
509  dp2.Assemble();
510  dp2.Finalize();
511 
512  HypreParMatrix * dP2 = dp2.ParallelAssemble();
513  temp = dP2->EliminateRowsCols(ess_tdof_list);
514  delete temp;
515 
516  Vector b2(y.Size());
517  dP2->Mult(y, b2);
518  delete dP2;
519 
520  double dp1_result = InnerProduct(pfes->GetComm(), yB, b1);
521  double dp2_result = InnerProduct(pfes->GetComm(), yB, b2);
522 
523  qBdot[0] = -dp1_result;
524  qBdot[1] = -dp2_result;
525 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:2052
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:551
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
Conjugate gradient method.
Definition: solvers.hpp:316
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:267
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:476
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
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:493
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:1931
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:150
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:9143
Parallel smoothers in hypre.
Definition: hypre.hpp:840
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
int SpaceDimension() const
Definition: mesh.hpp:912
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:99
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:204
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:702
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:316
void InitB(TimeDependentAdjointOperator &f_)
Initialize the adjoint problem.
Definition: sundials.cpp:862
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
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:828
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:1529
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:277
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:3028
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)