MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
cvsRoberts_ASAi_dns.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
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//
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
81using namespace std;
82using namespace mfem;
83
84
85// We create a TimeDependentAdjointOperator implementation of the rate equations
86// to recreate the cvsRoberts_ASAi_dns problem
88{
89public:
90 RobertsTDAOperator(int dim, Vector p) :
92 p_(p),
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 real_t t, const Vector &y,
113 const Vector &yB, const Vector &fyB, int jokB,
114 int *jcurB, real_t gammaB);
115
116 // Setup custom MFEM solve
117 virtual int SUNImplicitSolveB(Vector &x, const Vector &b, real_t tol);
118
119 ~RobertsTDAOperator()
120 {
122 }
123
124protected:
125 Vector p_;
126
127 // Solvers
130};
131
132int main(int argc, char *argv[])
133{
134 // Parse command-line options.
135 real_t t_final = 4e7;
136 real_t dt = 0.01;
137
138 // Relative tolerance for CVODES.
139 const real_t 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
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
173
174 // Set the initial time
175 real_t t = 0.0;
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)
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 real_t 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
217
218 // Initialize the adjoint solve
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 real_t 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;
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 real_t TBout1 = 40.;
251 Vector dG_dp(3);
252 dG_dp=0.;
253 if (cvodes)
254 {
255 t = t_final;
259 cvodes->SetMaxNStepsB(5000);
260
261 // Results at time TBout1
262 real_t 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
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
297void 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
306 Vector &qdot) const
307{
308 qdot[0] = y[2];
309}
310
312void RobertsTDAOperator::AdjointRateMult(const Vector &y, Vector & yB,
313 Vector &yBdot) const
314{
315 real_t l21 = (yB[1]-yB[0]);
316 real_t l32 = (yB[2]-yB[1]);
317 real_t p1 = p_[0];
318 real_t p2 = p_[1];
319 real_t 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
327 const Vector &yB,
328 Vector &qBdot) const
329{
330 real_t l21 = (yB[1]-yB[0]);
331 real_t l32 = (yB[2]-yB[1]);
332 real_t 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
340int RobertsTDAOperator::SUNImplicitSetupB(const real_t t, const Vector &y,
341 const Vector &yB,
342 const Vector &fyB, int jokB,
343 int *jcurB, real_t gammaB)
344{
345 // M = I- gamma J
346 // J = dfB/dyB
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.;
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;
366
367 return 0;
368}
369
370// cvsRoberts_ASAi_dns implicit solve for adjoint
371int RobertsTDAOperator::SUNImplicitSolveB(Vector &x, const Vector &b,
372 real_t tol)
373{
374 // The argument "tol" is ignored in this implementation for simplicity.
377 return (0);
378}
void SetMaxNStepsB(int mxstepsB)
Set the maximum number of backward steps.
void GetForwardSolution(double tB, mfem::Vector &yy)
Get Interpolated Forward solution y at backward integration time tB.
virtual void Step(Vector &x, double &t, double &dt)
void SetWFTolerances(EWTFunction func)
Set multiplicative error weights.
virtual void StepB(Vector &w, double &t, double &dt)
void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB=1e-3, double abstolQB=1e-8)
void InitQuadIntegration(mfem::Vector &q0, double reltolQ=1e-3, double abstolQ=1e-8)
Interface to the CVODE library – linear multi-step methods.
Definition sundials.hpp:386
void PrintInfo() const
Print various CVODE statistics.
Definition sundials.cpp:919
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition sundials.cpp:855
void SetSVtolerances(double reltol, Vector abstol)
Set the scalar relative and vector of absolute tolerances.
Definition sundials.cpp:881
GMRES method.
Definition solvers.hpp:547
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the GMRES method.
Definition solvers.cpp:980
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:179
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Data type sparse matrix.
Definition sparsemat.hpp:51
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR.
void Set(const int i, const int j, const real_t val)
real_t t
Current time.
Definition operator.hpp:340
virtual void SetTime(const real_t t_)
Set the current time.
Definition operator.hpp:360
Vector data type.
Definition vector.hpp:80
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition vector.cpp:755
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
int dim
Definition ex24.cpp:53
int main()
real_t b
Definition lissajous.cpp:42
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
real_t p(const Vector &x, real_t t)
RefCoord t[3]