MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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// Advection-Diffusion Miniapp: Parallel MFEM CVODES Example
14// ----------------------------------------------------------
15//
17//
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
77using namespace std;
78using namespace mfem;
79
82{
83public:
84 AdvDiffSUNDIALS(int ydot_dim, int ybdot_dim, Vector p,
85 ParFiniteElementSpace *fes, Array<int> & ess_tdof) :
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);
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);
113 k->Assemble(skip_zeros);
114 k->Finalize(skip_zeros);
115
116 k1 = new ParBilinearForm(pfes);
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
132 delete temp;
133
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 real_t gamma);
153
154 virtual int SUNImplicitSolve(const Vector &b, Vector &x, real_t tol);
155
156 virtual void QuadratureSensitivityMult(const Vector &y, const Vector &yB,
157 Vector &qbdot) const;
158
160 {
161 delete m;
162 delete k;
163 delete k1;
164 delete M;
165 delete K;
167 delete Mf;
168 delete p0;
169 delete mp0;
170 delete p2;
171 }
172
173protected:
174 Vector p_;
175 Array<int> ess_tdof_list;
177
178 // Internal matrices
181 ParBilinearForm *k1;
182
186
187 HypreParMatrix *Mf;
188
189 CGSolver M_solver;
190 HypreSmoother M_prec;
194};
195
196// Initial conditions for the problem
198{
199 return x[0]*(2. - x[0])*exp(2.*x[0]);
200}
201
202int 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 real_t t_final = 2.5;
213 real_t 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");
224 "Number of times to refine the mesh uniformly.");
226 "A switch to toggle between CV_ADAMS, and CV_BDF stepping modes");
228 "Final time; start time is 0.");
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
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
303
304 // Set the initial time to the TimeDependentAdjointOperator
305 real_t t = 0.0;
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);
312 cvodes->UseSundialsLinearSolver();
313 cvodes->SetMaxNSteps(5000);
314
315 // Relative and absolute tolerances for CVODES
316 real_t 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
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 real_t 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 real_t g = obj(u);
359 if (myid == 0)
360 {
361 cout << "g: " << g << endl;
362 }
363
365 ParGridFunction v(fes);
366 v = 1.;
367 v.SetSubVector(ess_tdof_list, 0.0);
369
370 // Initialize quadrature sensitivity values to zero
371 Vector qBdot(p.Size());
372 qBdot = 0.;
373
374 t = t_final;
377 cvodes->UseSundialsLinearSolverB();
378 cvodes->SetSStolerancesB(reltol, abstol);
379
380 // Results at time TBout1
381 real_t 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
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
414void 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
430 const Vector &fy,
431 int jok, int *jcur, real_t 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
445int AdvDiffSUNDIALS::SUNImplicitSolve(const Vector &b, Vector &x, real_t tol)
446{
447 Vector z(b.Size());
448 M->Mult(b,z);
449
450 CGSolver solver(pfes->GetComm());
451 HypreSmoother prec;
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
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);
472 M_solver.Mult(z, yBdot);
473}
474
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.);
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);
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 real_t dp1_result = InnerProduct(pfes->GetComm(), yB, b1);
518 real_t dp2_result = InnerProduct(pfes->GetComm(), yB, b2);
519
520 qBdot[0] = -dp1_result;
521 qBdot[1] = -dp2_result;
522}
real_t u_init(const Vector &x)
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void Print(std::ostream &out=mfem::out, int width=4) const
Prints array to stream with width elements per row.
Definition array.cpp:23
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
Adds new Domain Integrator. Assumes ownership of bfi.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
virtual void Step(Vector &x, double &t, double &dt)
void UseSundialsLinearSolverB()
Use built in SUNDIALS Newton solver.
virtual void StepB(Vector &w, double &t, double &dt)
void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB=1e-3, double abstolQB=1e-8)
void SetSStolerancesB(double reltol, double abstol)
Tolerance specification functions for the adjoint problem.
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition sundials.cpp:875
void SetMaxNSteps(int steps)
Set the maximum number of time steps.
Definition sundials.cpp:899
void PrintInfo() const
Print various CVODE statistics.
Definition sundials.cpp:919
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition sundials.cpp:855
A coefficient that is constant across space and time.
Class for domain integration .
Definition lininteg.hpp:109
A general function coefficient.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition hypre.cpp:2356
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, real_t alpha=1.0, real_t beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition hypre.cpp:1815
Wrapper for hypre's parallel vector class.
Definition hypre.hpp:206
void Print(const char *fname) const
Prints the locally owned rows in parallel.
Definition hypre.cpp:413
Parallel smoothers in hypre.
Definition hypre.hpp:1020
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition hypre.cpp:3550
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Adds new Domain Integrator. Assumes ownership of lfi.
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
static Mesh MakeCartesian1D(int n, real_t sx=1.0)
Creates 1D mesh, divided into n equal intervals.
Definition mesh.cpp:4235
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
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.
Class for parallel bilinear form.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:273
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
Class for parallel grid function.
Definition pgridfunc.hpp:33
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Class for parallel linear form.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Class for parallel meshes.
Definition pmesh.hpp:34
virtual void SetTime(const real_t t_)
Set the current time.
Definition operator.hpp:360
Vector coefficient that is constant in space and time.
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
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
real_t InnerProduct(HypreParVector *x, HypreParVector *y)
Definition hypre.cpp:439
float real_t
Definition config.hpp:43
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
real_t p(const Vector &x, real_t t)
RefCoord t[3]