MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
darcy.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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// Poisson/Darcy Mixed Method Solver
14// ---------------------------------
15//
16// Solves a Poisson problem -Delta p = f using a mixed finite element
17// formulation). The right-hand side of the Poisson problem is the same as that
18// used in the LOR Solvers miniapp (see miniapps/solvers). Dirichlet boundary
19// conditions are enforced on all domain boundaries.
20//
21// Optionally, the equation alpha*p - Delta p = f can be solved by setting the
22// alpha parameter to a nonzero value.
23
24// This can be written in the form of a Darcy problem
25//
26// -u - grad(p) = 0
27// alpha*p + div(u) = f
28//
29// where natural boundary conditions are enforced on the flux u, and the
30// Dirichlet condition on p is enforced by modifying the right-hand side.
31//
32// The resulting saddle-point system is solved using MINRES with a matrix-free
33// block-diagonal preconditioner.
34//
35// See also example 5 and its parallel version.
36//
37// Sample runs:
38//
39// darcy
40// mpirun -np 4 darcy -m ../../data/fichera-q2.mesh
41
42#include "mfem.hpp"
43#include <iostream>
44#include <memory>
45
48
50
51using namespace std;
52using namespace mfem;
53
54ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
55
56int main(int argc, char *argv[])
57{
58 Mpi::Init(argc, argv);
60
61 const char *mesh_file = "../../data/star.mesh";
62 const char *device_config = "cpu";
63 int ser_ref = 1;
64 int par_ref = 1;
65 int order = 3;
66 real_t alpha = 0.0;
67
68 OptionsParser args(argc, argv);
69 args.AddOption(&device_config, "-d", "--device",
70 "Device configuration string, see Device::Configure().");
71 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
72 args.AddOption(&ser_ref, "-rs", "--serial-refine",
73 "Number of times to refine the mesh in serial.");
74 args.AddOption(&par_ref, "-rp", "--parallel-refine",
75 "Number of times to refine the mesh in parallel.");
76 args.AddOption(&order, "-o", "--order", "Polynomial degree.");
77 args.AddOption(&alpha, "-a", "--alpha", "Value of alpha coefficient.");
78 args.ParseCheck();
79
80 Device device(device_config);
81 if (Mpi::Root()) { device.Print(); }
82
83 ParMesh mesh = LoadParMesh(mesh_file, ser_ref, par_ref);
84 const int dim = mesh.Dimension();
85 MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
86
88 const int mt = FiniteElement::VALUE;
89 RT_FECollection fec_rt(order-1, dim, b1, b2);
90 L2_FECollection fec_l2(order-1, dim, b2, mt);
91 ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
92 ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
93
94 HYPRE_BigInt ndofs_rt = fes_rt.GlobalTrueVSize();
95 HYPRE_BigInt ndofs_l2 = fes_l2.GlobalTrueVSize();
96
97 if (Mpi::Root())
98 {
99 cout << "\nRT DOFs: " << ndofs_rt << "\nL2 DOFs: " << ndofs_l2 << endl;
100 }
101
102 Array<int> ess_rt_dofs; // empty
103
104 // f is the RHS, u is the exact solution
105 FunctionCoefficient f_coeff(f(alpha)), u_coeff(u);
106
107 // Assemble the right-hand side for the scalar (L2) unknown.
108 ParLinearForm b_l2(&fes_l2);
109 b_l2.AddDomainIntegrator(new DomainLFIntegrator(f_coeff));
110 b_l2.UseFastAssembly(true);
111 b_l2.Assemble();
112
113 // Enforce Dirichlet boundary conditions on the scalar unknown by adding
114 // the boundary term to the flux equation.
115 ParLinearForm b_rt(&fes_rt);
117 b_rt.UseFastAssembly(true);
118 b_rt.Assemble();
119
120 if (Mpi::Root()) { cout << "\nSaddle point solver... " << flush; }
122
123 // Set up the block system of the form
124 //
125 // [ W D ][ u ] = [ f ]
126 // [ D^T -M ][ q ] = [ g_D ]
127 //
128 // where W is the L2 mass matrix, D is the discrete divergence, and M is
129 // the RT mass matrix.
130 //
131 // If the coefficient alpha is set to zero, the system takes the form
132 //
133 // [ 0 D ][ u ] = [ f ]
134 // [ D^T -M ][ q ] = [ g_D ]
135 //
136 // u is the scalar unknown, and q is the flux. f is the right-hand side from
137 // the Poisson problem, and g_D is the contribution to the right-hand side
138 // from the Dirichlet boundary condition.
139
140 ConstantCoefficient one(1.0);
141 ConstantCoefficient alpha_coeff(alpha);
142 const auto solver_mode = HdivSaddlePointSolver::Mode::DARCY;
143 HdivSaddlePointSolver saddle_point_solver(
144 mesh, fes_rt, fes_l2, alpha_coeff, one, ess_rt_dofs, solver_mode);
145
146 const Array<int> &offsets = saddle_point_solver.GetOffsets();
147 BlockVector X_block(offsets), B_block(offsets);
148
149 b_l2.ParallelAssemble(B_block.GetBlock(0));
150 b_rt.ParallelAssemble(B_block.GetBlock(1));
151 B_block.SyncFromBlocks();
152
153 X_block = 0.0;
154 saddle_point_solver.Mult(B_block, X_block);
155 X_block.SyncToBlocks();
156
157 if (Mpi::Root())
158 {
159 cout << "Done.\nIterations: "
160 << saddle_point_solver.GetNumIterations()
161 << "\nElapsed: " << tic_toc.RealTime() << endl;
162 }
163
164 ParGridFunction x(&fes_l2);
165 x.SetFromTrueDofs(X_block.GetBlock(0));
166 const real_t error = x.ComputeL2Error(u_coeff);
167 if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
168
169 return 0;
170}
171
172ParMesh LoadParMesh(const char *mesh_file, int ser_ref, int par_ref)
173{
174 Mesh serial_mesh = Mesh::LoadFromFile(mesh_file);
175 for (int i = 0; i < ser_ref; ++i) { serial_mesh.UniformRefinement(); }
176 ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
177 serial_mesh.Clear();
178 for (int i = 0; i < par_ref; ++i) { mesh.UniformRefinement(); }
179 return mesh;
180}
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
@ GaussLegendre
Open type.
Definition fe_base.hpp:31
A class to handle Vectors in a block fashion.
void SyncFromBlocks() const
Synchronize the memory location flags (i.e. the memory validity flags) of the big/monolithic block-ve...
Vector & GetBlock(int i)
Get the i-th vector in the block.
A coefficient that is constant across space and time.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
Class for domain integration .
Definition lininteg.hpp:109
A general function coefficient.
Solve the H(div) saddle-point system using MINRES with matrix-free block-diagonal preconditioning.
const Array< int > & GetOffsets() const
Return the offsets of the block system.
void Mult(const Vector &b, Vector &x) const override
Solve the linear system for L2 (scalar) and RT (flux) unknowns.
int GetNumIterations() const
Get the number of MINRES iterations.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator. Assumes ownership of lfi.
void UseFastAssembly(bool use_fa)
Which assembly algorithm to use: the new device-compatible fast assembly (true), or the legacy CPU-on...
Mesh data type.
Definition mesh.hpp:56
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
static Mesh LoadFromFile(const std::string &filename, int generate_edges=0, int refine=1, bool fix_orientation=true)
Definition mesh.cpp:4225
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void ParseCheck(std::ostream &out=mfem::out)
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
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
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void SetFromTrueDofs(const Vector &tv) override
Set the GridFunction from the given true-dof vector.
Class for parallel linear form.
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Class for parallel meshes.
Definition pmesh.hpp:34
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition tic_toc.cpp:432
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition tic_toc.cpp:411
void Clear()
Clear the elapsed time on the stopwatch and restart it if it's running.
Definition tic_toc.cpp:406
ParMesh LoadParMesh(const char *mesh_file, int ser_ref=0, int par_ref=0)
Definition darcy.cpp:172
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
StopWatch tic_toc
Definition tic_toc.cpp:450
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30