MFEM  v4.6.0
Finite element discretization library
cvsRoberts_ASAi_dns.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // cvsRoberts-ASAi-DNS Miniapp: Serial MFEM CVODES Example
14 // --------------------------------------------------------
15 //
16 // Compile with: make cvsRoberts_ASAi_dns
17 //
18 // Sample runs: cvsRoberts_ASAi_dns -dt 0.01
19 // cvsRoberts_ASAi_dns -dt 0.005
20 //
21 // Description: This example is a port of cvodes/serial/cvsRoberts_ASAi_dns
22 // example that is part of SUNDIALS. The goal is to demonstrate
23 // how to use the adjoint SUNDIALS CVODES interface from MFEM.
24 // Below is an excerpt description from the aforementioned file.
25 //
26 // Adjoint sensitivity example problem:
27 //
28 // The following is a simple example problem, with the coding needed for its
29 // solution by CVODES. The problem is from chemical kinetics, and consists of
30 // the following three rate equations
31 //
32 // dy1/dt = -p1*y1 + p2*y2*y3
33 // dy2/dt = p1*y1 - p2*y2*y3 - p3*(y2)^2
34 // dy3/dt = p3*(y2)^2
35 //
36 // on the interval from t = 0.0 to t = 4e10, with initial conditions: y1 = 1.0,
37 // y2 = y3 = 0. The reaction rates are: p1=0.04, p2=1e4, and p3=3e7. The problem
38 // is stiff. This program solves the problem with the BDF method, Newton
39 // iteration with the DENSE linear solver, and a user-supplied Jacobian routine.
40 // It uses a scalar relative tolerance and a vector absolute tolerance. Output
41 // is printed in decades from t = 0.4 to t = 4e10. Run statistics (optional
42 // outputs) are printed at the end.
43 //
44 // Optionally, CVODES can compute sensitivities with respect to the problem
45 // parameters p1, p2, and p3 of the following quantity:
46 //
47 // G = int_t0^t1 g(t,p,y) dt
48 //
49 // where g(t,p,y) = y3.
50 //
51 // The gradient dG/dp is obtained as:
52 //
53 // dG/dp = int_t0^t1 (g_p - lambda^T f_p ) dt - lambda^T(t0)*y0_p
54 // = - xi^T(t0) - lambda^T(t0)*y0_p
55 //
56 // where lambda and xi are solutions of:
57 //
58 // d(lambda)/dt = - (f_y)^T * lambda - (g_y)^T
59 // lambda(t1) = 0
60 //
61 // and
62 //
63 // d(xi)/dt = - (f_p)^T * lambda + (g_p)^T
64 // xi(t1) = 0
65 //
66 // During the backward integration, CVODES also evaluates G as
67 //
68 // G = - phi(t0)
69 //
70 // where d(phi)/dt = g(t,y,p), phi(t1) = 0.
71 
72 #include "mfem.hpp"
73 #include <fstream>
74 #include <iostream>
75 #include <algorithm>
76 
77 #ifndef MFEM_USE_SUNDIALS
78 #error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
79 #endif
80 
81 using namespace std;
82 using namespace mfem;
83 
84 
85 // We create a TimeDependentAdjointOperator implementation of the rate equations
86 // to recreate the cvsRoberts_ASAi_dns problem
87 class RobertsTDAOperator : public TimeDependentAdjointOperator
88 {
89 public:
90  RobertsTDAOperator(int dim, Vector p) :
92  p_(p),
93  adjointMatrix(NULL)
94  {}
95 
96  // Rate equation for forward problem
97  virtual void Mult(const Vector &x, Vector &y) const;
98 
99  // Quadrature integration for G
100  virtual void QuadratureIntegration(const Vector &x, Vector &y) const;
101 
102  // Adjoint rate equation corresponding to d(lambda)/dt
103  virtual void AdjointRateMult(const Vector &y, Vector &yB,
104  Vector &yBdot) const;
105 
106  // Quadrature sensitivity equations corresponding to dG/dp
107  virtual void QuadratureSensitivityMult(const Vector &y, const Vector &yB,
108  Vector &qbdot) const;
109 
110  // Setup custom MFEM solvers using GMRES since the Jacobian matrix is not
111  // symmetric
112  virtual int SUNImplicitSetupB(const double t, const Vector &y,
113  const Vector &yB, const Vector &fyB, int jokB,
114  int *jcurB, double gammaB);
115 
116  // Setup custom MFEM solve
117  virtual int SUNImplicitSolveB(Vector &x, const Vector &b, double tol);
118 
119  ~RobertsTDAOperator()
120  {
121  delete adjointMatrix;
122  }
123 
124 protected:
125  Vector p_;
126 
127  // Solvers
128  GMRESSolver adjointSolver;
129  SparseMatrix* adjointMatrix;
130 };
131 
132 int main(int argc, char *argv[])
133 {
134  // Parse command-line options.
135  double t_final = 4e7;
136  double dt = 0.01;
137 
138  // Relative tolerance for CVODES.
139  const double reltol = 1e-4;
140 
141  int precision = 8;
142  cout.precision(precision);
143 
144  OptionsParser args(argc, argv);
145  args.AddOption(&dt, "-dt", "--time-step", "Time step.");
146 
147  args.Parse();
148  if (!args.Good())
149  {
150  args.PrintUsage(cout);
151  return 1;
152  }
153  args.PrintOptions(cout);
154 
155  // The original cvsRoberts_ASAi_dns problem is a fixed sized problem of size
156  // 3. Define the solution vector.
157  Vector u(3);
158  u = 0.;
159  u[0] = 1.;
160 
161  // Define the TimeDependentAdjointOperator which implements the various rate
162  // equations: rate, quadrature rate, adjoint rate, and quadrature sensitivity
163  // rate equation
164 
165  // Define material parameters p
166  Vector p(3);
167  p[0] = 0.04;
168  p[1] = 1.0e4;
169  p[2] = 3.0e7;
170 
171  // 3 is the size of the adjoint solution vector
172  RobertsTDAOperator adv(3, p);
173 
174  // Set the initial time
175  double t = 0.0;
176  adv.SetTime(t);
177 
178  // Create the CVODES solver and set the various tolerances. Set absolute
179  // tolerances for the solution
180  Vector abstol_v(3);
181  abstol_v[0] = 1.0e-8;
182  abstol_v[1] = 1.0e-14;
183  abstol_v[2] = 1.0e-6;
184 
185  // Initialize the quadrature result
186  Vector q(1);
187  q = 0.;
188 
189  // Create the solver
190  CVODESSolver *cvodes = new CVODESSolver(CV_BDF);
191 
192  // Initialize the forward problem (this must be done first)
193  cvodes->Init(adv);
194 
195  // Set error control function
196  cvodes->SetWFTolerances([reltol, abstol_v]
197  (Vector y, Vector w, CVODESolver * self)
198  {
199  for (int i = 0; i < y.Size(); i++)
200  {
201  double ww = reltol * abs(y[i]) + abstol_v[i];
202  if (ww <= 0.) { return -1; }
203  w[i] = 1./ww;
204  }
205  return 0;
206  }
207  );
208 
209  // Set weighted tolerances
210  cvodes->SetSVtolerances(reltol, abstol_v);
211 
212  // Use the builtin sundials solver for the forward solver
213  cvodes->UseSundialsLinearSolver();
214 
215  // Initialize the quadrature integration and set the tolerances
216  cvodes->InitQuadIntegration(q, 1.e-6, 1.e-6);
217 
218  // Initialize the adjoint solve
219  cvodes->InitAdjointSolve(150, CV_HERMITE);
220 
221  // Perform time-integration (looping over the time iterations, ti, with a
222  // time-step dt).
223  bool done = false;
224  for (int ti = 0; !done; )
225  {
226  double dt_real = max(dt, t_final - t);
227  cvodes->Step(u, t, dt_real);
228  ti++;
229 
230  done = (t >= t_final - 1e-8*dt);
231 
232  if (done)
233  {
234  cvodes->PrintInfo();
235  }
236  }
237 
238  cout << "Final Solution: " << t << endl;
239  u.Print();
240 
241  q = 0.;
242  cout << " Final Quadrature " << endl;
243  cvodes->EvalQuadIntegration(t, q);
244  q.Print();
245 
246  // Solve the adjoint problem at different points in time. Create the adjoint
247  // solution vector
248  Vector w(3);
249  w=0.;
250  double TBout1 = 40.;
251  Vector dG_dp(3);
252  dG_dp=0.;
253  if (cvodes)
254  {
255  t = t_final;
256  adv.SetTime(t);
257  cvodes->InitB(adv);
258  cvodes->InitQuadIntegrationB(dG_dp, 1.e-6, 1.e-6);
259  cvodes->SetMaxNStepsB(5000);
260 
261  // Results at time TBout1
262  double dt_real = max(dt, t - TBout1);
263  cvodes->StepB(w, t, dt_real);
264  cout << "t: " << t << endl;
265  cout << "w:" << endl;
266  w.Print();
267 
268  cvodes->GetForwardSolution(t, u);
269  cout << "u:" << endl;
270  u.Print();
271 
272  // Results at T0
273  dt_real = max(dt, t - 0.);
274  cvodes->StepB(w, t, dt_real);
275  cout << "t: " << t << endl;
276  cout << "w:" << endl;
277 
278  cvodes->GetForwardSolution(t, u);
279  w.Print();
280  cout << "u:" << endl;
281  u.Print();
282 
283  // Evaluate Sensitivity
284  cvodes->EvalQuadIntegrationB(t, dG_dp);
285  cout << "dG/dp:" << endl;
286  dG_dp.Print();
287 
288  }
289 
290  // Free the used memory.
291  delete cvodes;
292 
293  return 0;
294 }
295 
296 // cvsRoberts_ASAi_dns rate equation
297 void RobertsTDAOperator::Mult(const Vector &x, Vector &y) const
298 {
299  y[0] = -p_[0]*x[0] + p_[1]*x[1]*x[2];
300  y[2] = p_[2]*x[1]*x[1];
301  y[1] = -y[0] - y[2];
302 }
303 
304 // cvsRoberts_ASAi_dns quadrature rate equation
305 void RobertsTDAOperator::QuadratureIntegration(const Vector &y,
306  Vector &qdot) const
307 {
308  qdot[0] = y[2];
309 }
310 
311 // cvsRoberts_ASAi_dns adjoint rate equation
312 void RobertsTDAOperator::AdjointRateMult(const Vector &y, Vector & yB,
313  Vector &yBdot) const
314 {
315  double l21 = (yB[1]-yB[0]);
316  double l32 = (yB[2]-yB[1]);
317  double p1 = p_[0];
318  double p2 = p_[1];
319  double p3 = p_[2];
320  yBdot[0] = -p1 * l21;
321  yBdot[1] = p2 * y[2] * l21 - 2. * p3 * y[1] * l32;
322  yBdot[2] = p2 * y[1] * l21 - 1.0;
323 }
324 
325 // cvsRoberts_ASAi_dns quadrature sensitivity rate equation
326 void RobertsTDAOperator::QuadratureSensitivityMult(const Vector &y,
327  const Vector &yB,
328  Vector &qBdot) const
329 {
330  double l21 = (yB[1]-yB[0]);
331  double l32 = (yB[2]-yB[1]);
332  double y23 = y[1] * y[2];
333 
334  qBdot[0] = y[0] * l21;
335  qBdot[1] = -y23 * l21;
336  qBdot[2] = y[1]*y[1]*l32;
337 }
338 
339 // cvsRoberts_ASAi_dns implicit solve setup for adjoint
340 int RobertsTDAOperator::SUNImplicitSetupB(const double t, const Vector &y,
341  const Vector &yB,
342  const Vector &fyB, int jokB,
343  int *jcurB, double gammaB)
344 {
345  // M = I- gamma J
346  // J = dfB/dyB
347  delete adjointMatrix;
348  adjointMatrix = new SparseMatrix(y.Size(), yB.Size());
349  for (int j = 0; j < y.Size(); j++)
350  {
351  Vector JacBj(yB.Size());
352  Vector yBone(yB.Size());
353  yBone = 0.;
354  yBone[j] = 1.;
355  AdjointRateMult(y, yBone, JacBj);
356  JacBj[2] += 1.;
357  for (int i = 0; i < y.Size(); i++)
358  {
359  adjointMatrix->Set(i,j, (i == j ? 1.0 : 0.) - gammaB * JacBj[i]);
360  }
361  }
362 
363  *jcurB = 1;
364  adjointMatrix->Finalize();
365  adjointSolver.SetOperator(*adjointMatrix);
366 
367  return 0;
368 }
369 
370 // cvsRoberts_ASAi_dns implicit solve for adjoint
371 int RobertsTDAOperator::SUNImplicitSolveB(Vector &x, const Vector &b,
372  double tol)
373 {
374  // The argument "tol" is ignored in this implementation for simplicity.
375  adjointSolver.SetRelTol(1e-14);
376  adjointSolver.Mult(b, x);
377  return (0);
378 }
void GetForwardSolution(double tB, mfem::Vector &yy)
Get Interpolated Forward solution y at backward integration time tB.
Definition: sundials.cpp:1025
virtual void StepB(Vector &w, double &t, double &dt)
Solve one adjoint time step.
Definition: sundials.cpp:1301
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:756
void SetSVtolerances(double reltol, Vector abstol)
Set the scalar relative and vector of absolute tolerances.
Definition: sundials.cpp:881
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
void EvalQuadIntegration(double t, Vector &q)
Evaluate Quadrature.
Definition: sundials.cpp:1005
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
virtual void Step(Vector &x, double &t, double &dt)
Definition: sundials.cpp:1273
void InitAdjointSolve(int steps, int interpolation)
Initialize Adjoint.
Definition: sundials.cpp:1072
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:50
Interface to the CVODE library – linear multi-step methods.
Definition: sundials.hpp:385
void InitQuadIntegration(mfem::Vector &q0, double reltolQ=1e-3, double abstolQ=1e-8)
Definition: sundials.cpp:1084
void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB=1e-3, double abstolQB=1e-8)
Initialize Quadrature Integration (Adjoint)
Definition: sundials.cpp:1099
double b
Definition: lissajous.cpp:42
void PrintInfo() const
Print various CVODE statistics.
Definition: sundials.cpp:919
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
GMRES method.
Definition: solvers.hpp:525
int dim
Definition: ex24.cpp:53
void InitB(TimeDependentAdjointOperator &f_)
Initialize the adjoint problem.
Definition: sundials.cpp:1039
void SetWFTolerances(EWTFunction func)
Set multiplicative error weights.
Definition: sundials.cpp:1211
int main(int argc, char *argv[])
RefCoord t[3]
Vector data type.
Definition: vector.hpp:58
void Init(TimeDependentAdjointOperator &f_)
Definition: sundials.cpp:1034
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
void SetMaxNStepsB(int mxstepsB)
Set the maximum number of backward steps.
Definition: sundials.cpp:1078
void EvalQuadIntegrationB(double t, Vector &dG_dp)
Evaluate Quadrature solution.
Definition: sundials.cpp:1015
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition: sundials.cpp:855