MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ex16p.cpp
Go to the documentation of this file.
1// MFEM Example 16 - Parallel Version
2//
3// Compile with: make ex16p
4//
5// Sample runs: mpirun -np 4 ex16p
6// mpirun -np 4 ex16p -m ../data/inline-tri.mesh
7// mpirun -np 4 ex16p -m ../data/disc-nurbs.mesh -tf 2
8// mpirun -np 4 ex16p -s 1 -a 0.0 -k 1.0
9// mpirun -np 4 ex16p -s 2 -a 1.0 -k 0.0
10// mpirun -np 8 ex16p -s 3 -a 0.5 -k 0.5 -o 4
11// mpirun -np 4 ex16p -s 14 -dt 1.0e-4 -tf 4.0e-2 -vs 40
12// mpirun -np 16 ex16p -m ../data/fichera-q2.mesh
13// mpirun -np 16 ex16p -m ../data/fichera-mixed.mesh
14// mpirun -np 16 ex16p -m ../data/escher-p2.mesh
15// mpirun -np 8 ex16p -m ../data/beam-tet.mesh -tf 10 -dt 0.1
16// mpirun -np 4 ex16p -m ../data/amr-quad.mesh -o 4 -rs 0 -rp 0
17// mpirun -np 4 ex16p -m ../data/amr-hex.mesh -o 2 -rs 0 -rp 0
18//
19// Description: This example solves a time dependent nonlinear heat equation
20// problem of the form du/dt = C(u), with a non-linear diffusion
21// operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u.
22//
23// The example demonstrates the use of nonlinear operators (the
24// class ConductionOperator defining C(u)), as well as their
25// implicit time integration. Note that implementing the method
26// ConductionOperator::ImplicitSolve is the only requirement for
27// high-order implicit (SDIRK) time integration. In this example,
28// the diffusion operator is linearized by evaluating with the
29// lagged solution from the previous timestep, so there is only
30// a linear solve. Optional saving with ADIOS2
31// (adios2.readthedocs.io) is also illustrated.
32//
33// We recommend viewing examples 2, 9 and 10 before viewing this
34// example.
35
36#include "mfem.hpp"
37#include <fstream>
38#include <iostream>
39
40using namespace std;
41using namespace mfem;
42
43/** After spatial discretization, the conduction model can be written as:
44 *
45 * du/dt = M^{-1}(-Ku)
46 *
47 * where u is the vector representing the temperature, M is the mass matrix,
48 * and K is the diffusion operator with diffusivity depending on u:
49 * (\kappa + \alpha u).
50 *
51 * Class ConductionOperator represents the right-hand side of the above ODE.
52 */
53class ConductionOperator : public TimeDependentOperator
54{
55protected:
56 ParFiniteElementSpace &fespace;
57 Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
58
61
62 HypreParMatrix Mmat;
63 HypreParMatrix Kmat;
64 HypreParMatrix *T; // T = M + dt K
65 real_t current_dt;
66
67 CGSolver M_solver; // Krylov solver for inverting the mass matrix M
68 HypreSmoother M_prec; // Preconditioner for the mass matrix M
69
70 CGSolver T_solver; // Implicit solver for T = M + dt K
71 HypreSmoother T_prec; // Preconditioner for the implicit solver
72
73 real_t alpha, kappa;
74
75 mutable Vector z; // auxiliary vector
76
77public:
78 ConductionOperator(ParFiniteElementSpace &f, real_t alpha, real_t kappa,
79 const Vector &u);
80
81 virtual void Mult(const Vector &u, Vector &du_dt) const;
82 /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
83 This is the only requirement for high-order SDIRK implicit integration.*/
84 virtual void ImplicitSolve(const real_t dt, const Vector &u, Vector &k);
85
86 /// Update the diffusion BilinearForm K using the given true-dof vector `u`.
87 void SetParameters(const Vector &u);
88
89 virtual ~ConductionOperator();
90};
91
93
94int main(int argc, char *argv[])
95{
96 // 1. Initialize MPI and HYPRE.
97 Mpi::Init(argc, argv);
98 int num_procs = Mpi::WorldSize();
99 int myid = Mpi::WorldRank();
100 Hypre::Init();
101
102 // 2. Parse command-line options.
103 const char *mesh_file = "../data/star.mesh";
104 int ser_ref_levels = 2;
105 int par_ref_levels = 1;
106 int order = 2;
107 int ode_solver_type = 3;
108 real_t t_final = 0.5;
109 real_t dt = 1.0e-2;
110 real_t alpha = 1.0e-2;
111 real_t kappa = 0.5;
112 bool visualization = true;
113 bool visit = false;
114 int vis_steps = 5;
115 bool adios2 = false;
116
117 int precision = 8;
118 cout.precision(precision);
119
120 OptionsParser args(argc, argv);
121 args.AddOption(&mesh_file, "-m", "--mesh",
122 "Mesh file to use.");
123 args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
124 "Number of times to refine the mesh uniformly in serial.");
125 args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
126 "Number of times to refine the mesh uniformly in parallel.");
127 args.AddOption(&order, "-o", "--order",
128 "Order (degree) of the finite elements.");
129 args.AddOption(&ode_solver_type, "-s", "--ode-solver",
130 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
131 "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
132 args.AddOption(&t_final, "-tf", "--t-final",
133 "Final time; start time is 0.");
134 args.AddOption(&dt, "-dt", "--time-step",
135 "Time step.");
136 args.AddOption(&alpha, "-a", "--alpha",
137 "Alpha coefficient.");
138 args.AddOption(&kappa, "-k", "--kappa",
139 "Kappa coefficient offset.");
140 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
141 "--no-visualization",
142 "Enable or disable GLVis visualization.");
143 args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
144 "--no-visit-datafiles",
145 "Save data files for VisIt (visit.llnl.gov) visualization.");
146 args.AddOption(&vis_steps, "-vs", "--visualization-steps",
147 "Visualize every n-th timestep.");
148 args.AddOption(&adios2, "-adios2", "--adios2-streams", "-no-adios2",
149 "--no-adios2-streams",
150 "Save data using adios2 streams.");
151 args.Parse();
152 if (!args.Good())
153 {
154 args.PrintUsage(cout);
155 return 1;
156 }
157
158 if (myid == 0)
159 {
160 args.PrintOptions(cout);
161 }
162
163 // 3. Read the serial mesh from the given mesh file on all processors. We can
164 // handle triangular, quadrilateral, tetrahedral and hexahedral meshes
165 // with the same code.
166 Mesh *mesh = new Mesh(mesh_file, 1, 1);
167 int dim = mesh->Dimension();
168
169 // 4. Define the ODE solver used for time integration. Several implicit
170 // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
171 // explicit Runge-Kutta methods are available.
172 ODESolver *ode_solver;
173 switch (ode_solver_type)
174 {
175 // Implicit L-stable methods
176 case 1: ode_solver = new BackwardEulerSolver; break;
177 case 2: ode_solver = new SDIRK23Solver(2); break;
178 case 3: ode_solver = new SDIRK33Solver; break;
179 // Explicit methods
180 case 11: ode_solver = new ForwardEulerSolver; break;
181 case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
182 case 13: ode_solver = new RK3SSPSolver; break;
183 case 14: ode_solver = new RK4Solver; break;
184 case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
185 // Implicit A-stable methods (not L-stable)
186 case 22: ode_solver = new ImplicitMidpointSolver; break;
187 case 23: ode_solver = new SDIRK23Solver; break;
188 case 24: ode_solver = new SDIRK34Solver; break;
189 default:
190 cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
191 delete mesh;
192 return 3;
193 }
194
195 // 5. Refine the mesh in serial to increase the resolution. In this example
196 // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
197 // a command-line parameter.
198 for (int lev = 0; lev < ser_ref_levels; lev++)
199 {
200 mesh->UniformRefinement();
201 }
202
203 // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
204 // this mesh further in parallel to increase the resolution. Once the
205 // parallel mesh is defined, the serial mesh can be deleted.
206 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
207 delete mesh;
208 for (int lev = 0; lev < par_ref_levels; lev++)
209 {
210 pmesh->UniformRefinement();
211 }
212
213 // 7. Define the vector finite element space representing the current and the
214 // initial temperature, u_ref.
215 H1_FECollection fe_coll(order, dim);
216 ParFiniteElementSpace fespace(pmesh, &fe_coll);
217
218 HYPRE_BigInt fe_size = fespace.GlobalTrueVSize();
219 if (myid == 0)
220 {
221 cout << "Number of temperature unknowns: " << fe_size << endl;
222 }
223
224 ParGridFunction u_gf(&fespace);
225
226 // 8. Set the initial conditions for u. All boundaries are considered
227 // natural.
229 u_gf.ProjectCoefficient(u_0);
230 Vector u;
231 u_gf.GetTrueDofs(u);
232
233 // 9. Initialize the conduction operator and the VisIt visualization.
234 ConductionOperator oper(fespace, alpha, kappa, u);
235
236 u_gf.SetFromTrueDofs(u);
237 {
238 ostringstream mesh_name, sol_name;
239 mesh_name << "ex16-mesh." << setfill('0') << setw(6) << myid;
240 sol_name << "ex16-init." << setfill('0') << setw(6) << myid;
241 ofstream omesh(mesh_name.str().c_str());
242 omesh.precision(precision);
243 pmesh->Print(omesh);
244 ofstream osol(sol_name.str().c_str());
245 osol.precision(precision);
246 u_gf.Save(osol);
247 }
248
249 VisItDataCollection visit_dc("Example16-Parallel", pmesh);
250 visit_dc.RegisterField("temperature", &u_gf);
251 if (visit)
252 {
253 visit_dc.SetCycle(0);
254 visit_dc.SetTime(0.0);
255 visit_dc.Save();
256 }
257
258 // Optionally output a BP (binary pack) file using ADIOS2. This can be
259 // visualized with the ParaView VTX reader.
260#ifdef MFEM_USE_ADIOS2
261 ADIOS2DataCollection* adios2_dc = NULL;
262 if (adios2)
263 {
264 std::string postfix(mesh_file);
265 postfix.erase(0, std::string("../data/").size() );
266 postfix += "_o" + std::to_string(order);
267 postfix += "_solver" + std::to_string(ode_solver_type);
268 const std::string collection_name = "ex16-p-" + postfix + ".bp";
269
270 adios2_dc = new ADIOS2DataCollection(MPI_COMM_WORLD, collection_name, pmesh);
271 adios2_dc->SetParameter("SubStreams", std::to_string(num_procs/2) );
272 adios2_dc->RegisterField("temperature", &u_gf);
273 adios2_dc->SetCycle(0);
274 adios2_dc->SetTime(0.0);
275 adios2_dc->Save();
276 }
277#endif
278
279 socketstream sout;
280 if (visualization)
281 {
282 char vishost[] = "localhost";
283 int visport = 19916;
284 sout.open(vishost, visport);
285 sout << "parallel " << num_procs << " " << myid << endl;
286 int good = sout.good(), all_good;
287 MPI_Allreduce(&good, &all_good, 1, MPI_INT, MPI_MIN, pmesh->GetComm());
288 if (!all_good)
289 {
290 sout.close();
291 visualization = false;
292 if (myid == 0)
293 {
294 cout << "Unable to connect to GLVis server at "
295 << vishost << ':' << visport << endl;
296 cout << "GLVis visualization disabled.\n";
297 }
298 }
299 else
300 {
301 sout.precision(precision);
302 sout << "solution\n" << *pmesh << u_gf;
303 sout << "pause\n";
304 sout << flush;
305 if (myid == 0)
306 {
307 cout << "GLVis visualization paused."
308 << " Press space (in the GLVis window) to resume it.\n";
309 }
310 }
311 }
312
313 // 10. Perform time-integration (looping over the time iterations, ti, with a
314 // time-step dt).
315 ode_solver->Init(oper);
316 real_t t = 0.0;
317
318 bool last_step = false;
319 for (int ti = 1; !last_step; ti++)
320 {
321 if (t + dt >= t_final - dt/2)
322 {
323 last_step = true;
324 }
325
326 ode_solver->Step(u, t, dt);
327
328 if (last_step || (ti % vis_steps) == 0)
329 {
330 if (myid == 0)
331 {
332 cout << "step " << ti << ", t = " << t << endl;
333 }
334
335 u_gf.SetFromTrueDofs(u);
336 if (visualization)
337 {
338 sout << "parallel " << num_procs << " " << myid << "\n";
339 sout << "solution\n" << *pmesh << u_gf << flush;
340 }
341
342 if (visit)
343 {
344 visit_dc.SetCycle(ti);
345 visit_dc.SetTime(t);
346 visit_dc.Save();
347 }
348
349#ifdef MFEM_USE_ADIOS2
350 if (adios2)
351 {
352 adios2_dc->SetCycle(ti);
353 adios2_dc->SetTime(t);
354 adios2_dc->Save();
355 }
356#endif
357 }
358 oper.SetParameters(u);
359 }
360
361#ifdef MFEM_USE_ADIOS2
362 if (adios2)
363 {
364 delete adios2_dc;
365 }
366#endif
367
368 // 11. Save the final solution in parallel. This output can be viewed later
369 // using GLVis: "glvis -np <np> -m ex16-mesh -g ex16-final".
370 {
371 ostringstream sol_name;
372 sol_name << "ex16-final." << setfill('0') << setw(6) << myid;
373 ofstream osol(sol_name.str().c_str());
374 osol.precision(precision);
375 u_gf.Save(osol);
376 }
377
378 // 12. Free the used memory.
379 delete ode_solver;
380 delete pmesh;
381
382 return 0;
383}
384
385ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, real_t al,
386 real_t kap, const Vector &u)
387 : TimeDependentOperator(f.GetTrueVSize(), (real_t) 0.0), fespace(f),
388 M(NULL), K(NULL), T(NULL), current_dt(0.0),
389 M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
390{
391 const real_t rel_tol = 1e-8;
392
393 M = new ParBilinearForm(&fespace);
394 M->AddDomainIntegrator(new MassIntegrator());
395 M->Assemble(0); // keep sparsity pattern of M and K the same
396 M->FormSystemMatrix(ess_tdof_list, Mmat);
397
398 M_solver.iterative_mode = false;
399 M_solver.SetRelTol(rel_tol);
400 M_solver.SetAbsTol(0.0);
401 M_solver.SetMaxIter(100);
402 M_solver.SetPrintLevel(0);
403 M_prec.SetType(HypreSmoother::Jacobi);
404 M_solver.SetPreconditioner(M_prec);
405 M_solver.SetOperator(Mmat);
406
407 alpha = al;
408 kappa = kap;
409
410 T_solver.iterative_mode = false;
411 T_solver.SetRelTol(rel_tol);
412 T_solver.SetAbsTol(0.0);
413 T_solver.SetMaxIter(100);
414 T_solver.SetPrintLevel(0);
415 T_solver.SetPreconditioner(T_prec);
416
417 SetParameters(u);
418}
419
420void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
421{
422 // Compute:
423 // du_dt = M^{-1}*-Ku
424 // for du_dt, where K is linearized by using u from the previous timestep
425 Kmat.Mult(u, z);
426 z.Neg(); // z = -z
427 M_solver.Mult(z, du_dt);
428}
429
430void ConductionOperator::ImplicitSolve(const real_t dt,
431 const Vector &u, Vector &du_dt)
432{
433 // Solve the equation:
434 // du_dt = M^{-1}*[-K(u + dt*du_dt)]
435 // for du_dt, where K is linearized by using u from the previous timestep
436 if (!T)
437 {
438 T = Add(1.0, Mmat, dt, Kmat);
439 current_dt = dt;
440 T_solver.SetOperator(*T);
441 }
442 MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
443 Kmat.Mult(u, z);
444 z.Neg();
445 T_solver.Mult(z, du_dt);
446}
447
448void ConductionOperator::SetParameters(const Vector &u)
449{
450 ParGridFunction u_alpha_gf(&fespace);
451 u_alpha_gf.SetFromTrueDofs(u);
452 for (int i = 0; i < u_alpha_gf.Size(); i++)
453 {
454 u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
455 }
456
457 delete K;
458 K = new ParBilinearForm(&fespace);
459
460 GridFunctionCoefficient u_coeff(&u_alpha_gf);
461
463 K->Assemble(0); // keep sparsity pattern of M and K the same
464 K->FormSystemMatrix(ess_tdof_list, Kmat);
465 delete T;
466 T = NULL; // re-compute T on the next ImplicitSolve
467}
468
469ConductionOperator::~ConductionOperator()
470{
471 delete T;
472 delete M;
473 delete K;
474}
475
477{
478 if (x.Norml2() < 0.5)
479 {
480 return 2.0;
481 }
482 else
483 {
484 return 1.0;
485 }
486}
void SetParameter(const std::string key, const std::string value) noexcept
Backward Euler ODE solver. L-stable.
Definition ode.hpp:412
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Conjugate gradient method.
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 RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
The classical forward Euler method.
Definition ode.hpp:118
A general function coefficient.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Parallel smoothers in hypre.
Definition hypre.hpp:1020
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Implicit midpoint method. A-stable, not L-stable.
Definition ode.hpp:425
Mesh data type.
Definition mesh.hpp:56
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
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 int WorldSize()
Return the size of 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).
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition ode.hpp:24
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition ode.cpp:18
virtual void Step(Vector &x, real_t &t, real_t &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
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.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
void SetFromTrueDofs(const Vector &tv) override
Set the GridFunction from the given true-dof vector.
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition ode.hpp:151
The classical explicit forth-order Runge-Kutta method, RK4.
Definition ode.hpp:164
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Base abstract class for first order time dependent operators.
Definition operator.hpp:317
Vector data type.
Definition vector.hpp:80
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
Data collection with VisIt I/O routines.
virtual void Save() override
Save the collection and a VisIt root file.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int close()
Close the socketstream.
const real_t alpha
Definition ex15.cpp:369
real_t InitialTemperature(const Vector &x)
Definition ex16p.cpp:476
real_t kappa
Definition ex24.cpp:54
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
const char vishost[]
RefCoord t[3]