MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
multidomain.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// This miniapp aims to demonstrate how to solve two PDEs, that represent
13// different physics, on the same domain. MFEM's SubMesh interface is used to
14// compute on and transfer between the spaces of predefined parts of the domain.
15// For the sake of simplicity, the spaces on each domain are using the same
16// order H1 finite elements. This does not mean that the approach is limited to
17// this configuration.
18//
19// A 3D domain comprised of an outer box with a cylinder shaped inside is used.
20//
21// A heat equation is described on the outer box domain
22//
23// dT/dt = κΔT in outer box
24// T = T_wall on outside wall
25// ∇T•n = 0 on inside (cylinder) wall
26//
27// with temperature T and coefficient κ (non-physical in this example).
28//
29// A convection-diffusion equation is described inside the cylinder domain
30//
31// dT/dt = κΔT - α∇•(b T) in inner cylinder
32// T = T_wall on cylinder wall (obtained from heat equation)
33// ∇T•n = 0 else
34//
35// with temperature T, coefficients κ, α and prescribed velocity profile b.
36//
37// To couple the solutions of both equations, a segregated solve with one way
38// coupling approach is used. The heat equation of the outer box is solved from
39// the timestep T_box(t) to T_box(t+dt). Then for the convection-diffusion
40// equation T_wall is set to T_box(t+dt) and the equation is solved for T(t+dt)
41// which results in a first-order one way coupling.
42
43#include "mfem.hpp"
44#include <fstream>
45#include <iostream>
46
47using namespace mfem;
48
49// Prescribed velocity profile for the convection-diffusion equation inside the
50// cylinder. The profile is constructed s.t. it approximates a no-slip (v=0)
51// directly at the cylinder wall boundary.
52void velocity_profile(const Vector &c, Vector &q)
53{
54 real_t A = 1.0;
55 real_t x = c(0);
56 real_t y = c(1);
57 real_t r = sqrt(pow(x, 2.0) + pow(y, 2.0));
58
59 q(0) = 0.0;
60 q(1) = 0.0;
61
62 if (std::abs(r) >= 0.25 - 1e-8)
63 {
64 q(2) = 0.0;
65 }
66 else
67 {
68 q(2) = A * exp(-(pow(x, 2.0) / 2.0 + pow(y, 2.0) / 2.0));
69 }
70}
71
72/**
73 * @brief Convection-diffusion time dependent operator
74 *
75 * dT/dt = κΔT - α∇•(b T)
76 *
77 * Can also be used to create a diffusion or convection only operator by setting
78 * α or κ to zero.
79 */
80class ConvectionDiffusionTDO : public TimeDependentOperator
81{
82public:
83 /**
84 * @brief Construct a new convection-diffusion time dependent operator.
85 *
86 * @param fes The ParFiniteElementSpace the solution is defined on
87 * @param ess_tdofs All essential true dofs (relevant if fes is using H1
88 * finite elements)
89 * @param alpha The convection coefficient
90 * @param kappa The diffusion coefficient
91 */
92 ConvectionDiffusionTDO(ParFiniteElementSpace &fes,
93 Array<int> ess_tdofs,
94 real_t alpha = 1.0,
95 real_t kappa = 1.0e-1)
96 : TimeDependentOperator(fes.GetTrueVSize()),
97 Mform(&fes),
98 Kform(&fes),
99 bform(&fes),
100 ess_tdofs_(ess_tdofs),
101 M_solver(fes.GetComm())
102 {
103 d = new ConstantCoefficient(-kappa);
106
108 Mform.Assemble(0);
109 Mform.Finalize();
110
111 if (fes.IsDGSpace())
112 {
113 M.Reset(Mform.ParallelAssemble(), true);
114
115 inflow = new ConstantCoefficient(0.0);
117 new BoundaryFlowIntegrator(*inflow, *q, alpha));
118 }
119 else
120 {
123 Kform.Assemble(0);
124
125 Array<int> empty;
126 Kform.FormSystemMatrix(empty, K);
127 Mform.FormSystemMatrix(ess_tdofs_, M);
128
129 bform.Assemble();
130 b = bform.ParallelAssemble();
131 }
132
133 M_solver.iterative_mode = false;
134 M_solver.SetRelTol(1e-8);
135 M_solver.SetAbsTol(0.0);
136 M_solver.SetMaxIter(100);
137 M_solver.SetPrintLevel(0);
139 M_solver.SetPreconditioner(M_prec);
140 M_solver.SetOperator(*M);
141
142 t1.SetSize(height);
143 t2.SetSize(height);
144 }
145
146 void Mult(const Vector &u, Vector &du_dt) const override
147 {
148 K->Mult(u, t1);
150 M_solver.Mult(t1, du_dt);
151 du_dt.SetSubVector(ess_tdofs_, 0.0);
152 }
153
154 ~ConvectionDiffusionTDO()
155 {
156 delete q;
157 delete d;
158 delete b;
159 }
160
161 /// Mass form
162 ParBilinearForm Mform;
163
164 /// Stiffness form. Might include diffusion, convection or both.
165 ParBilinearForm Kform;
166
167 /// Mass opeperator
169
170 /// Stiffness opeperator. Might include diffusion, convection or both.
172
173 /// RHS form
174 ParLinearForm bform;
175
176 /// RHS vector
177 Vector *b = nullptr;
178
179 /// Velocity coefficient
180 VectorCoefficient *q = nullptr;
181
182 /// Diffusion coefficient
183 Coefficient *d = nullptr;
184
185 /// Inflow coefficient
186 Coefficient *inflow = nullptr;
187
188 /// Essential true dof array. Relevant for eliminating boundary conditions
189 /// when using an H1 space.
190 Array<int> ess_tdofs_;
191
192 real_t current_dt = -1.0;
193
194 /// Mass matrix solver
195 CGSolver M_solver;
196
197 /// Mass matrix preconditioner
198 HypreSmoother M_prec;
199
200 /// Auxiliary vectors
201 mutable Vector t1, t2;
202};
203
204int main(int argc, char *argv[])
205{
206 Mpi::Init();
207 Hypre::Init();
208 int num_procs = Mpi::WorldSize();
209 int myid = Mpi::WorldRank();
210
211 int order = 2;
212 real_t t_final = 5.0;
213 real_t dt = 1.0e-5;
214 bool visualization = true;
215 int vis_steps = 10;
216
217 OptionsParser args(argc, argv);
219 "Finite element order (polynomial degree).");
221 "Final time; start time is 0.");
223 "Time step.");
225 "--no-visualization",
226 "Enable or disable GLVis visualization.");
228 "Visualize every n-th timestep.");
229 args.ParseCheck();
230
231 Mesh *serial_mesh = new Mesh("multidomain-hex.mesh");
232 ParMesh parent_mesh = ParMesh(MPI_COMM_WORLD, *serial_mesh);
233 delete serial_mesh;
234
235 parent_mesh.UniformRefinement();
236
237 H1_FECollection fec(order, parent_mesh.Dimension());
238
239 // Create the sub-domains and accompanying Finite Element spaces from
240 // corresponding attributes. This specific mesh has two domain attributes and
241 // 9 boundary attributes.
242 Array<int> cylinder_domain_attributes(1);
243 cylinder_domain_attributes[0] = 1;
244
245 auto cylinder_submesh =
246 ParSubMesh::CreateFromDomain(parent_mesh, cylinder_domain_attributes);
247
248 ParFiniteElementSpace fes_cylinder(&cylinder_submesh, &fec);
249
250 Array<int> inflow_attributes(cylinder_submesh.bdr_attributes.Max());
251 inflow_attributes = 0;
252 inflow_attributes[7] = 1;
253
254 Array<int> inner_cylinder_wall_attributes(
255 cylinder_submesh.bdr_attributes.Max());
256 inner_cylinder_wall_attributes = 0;
257 inner_cylinder_wall_attributes[8] = 1;
258
259 // For the convection-diffusion equation inside the cylinder domain, the
260 // inflow surface and outer wall are treated as Dirichlet boundary
261 // conditions.
262 Array<int> inflow_tdofs, interface_tdofs, ess_tdofs;
263 fes_cylinder.GetEssentialTrueDofs(inflow_attributes, inflow_tdofs);
264 fes_cylinder.GetEssentialTrueDofs(inner_cylinder_wall_attributes,
265 interface_tdofs);
266 ess_tdofs.Append(inflow_tdofs);
267 ess_tdofs.Append(interface_tdofs);
268 ess_tdofs.Sort();
269 ess_tdofs.Unique();
270 ConvectionDiffusionTDO cd_tdo(fes_cylinder, ess_tdofs);
271
272 ParGridFunction temperature_cylinder_gf(&fes_cylinder);
273 temperature_cylinder_gf = 0.0;
274
275 Vector temperature_cylinder;
276 temperature_cylinder_gf.GetTrueDofs(temperature_cylinder);
277
278 RK3SSPSolver cd_ode_solver;
279 cd_ode_solver.Init(cd_tdo);
280
281 Array<int> outer_domain_attributes(1);
282 outer_domain_attributes[0] = 2;
283
284 auto block_submesh = ParSubMesh::CreateFromDomain(parent_mesh,
285 outer_domain_attributes);
286
287 ParFiniteElementSpace fes_block(&block_submesh, &fec);
288
289 Array<int> block_wall_attributes(block_submesh.bdr_attributes.Max());
290 block_wall_attributes = 0;
291 block_wall_attributes[0] = 1;
292 block_wall_attributes[1] = 1;
293 block_wall_attributes[2] = 1;
294 block_wall_attributes[3] = 1;
295
296 Array<int> outer_cylinder_wall_attributes(
297 block_submesh.bdr_attributes.Max());
298 outer_cylinder_wall_attributes = 0;
299 outer_cylinder_wall_attributes[8] = 1;
300
301 fes_block.GetEssentialTrueDofs(block_wall_attributes, ess_tdofs);
302
303 ConvectionDiffusionTDO d_tdo(fes_block, ess_tdofs, 0.0, 1.0);
304
305 ParGridFunction temperature_block_gf(&fes_block);
306 temperature_block_gf = 0.0;
307
308 ConstantCoefficient one(1.0);
309 temperature_block_gf.ProjectBdrCoefficient(one, block_wall_attributes);
310
311 Vector temperature_block;
312 temperature_block_gf.GetTrueDofs(temperature_block);
313
314 RK3SSPSolver d_ode_solver;
315 d_ode_solver.Init(d_tdo);
316
317 Array<int> cylinder_surface_attributes(1);
318 cylinder_surface_attributes[0] = 9;
319
320 auto cylinder_surface_submesh = ParSubMesh::CreateFromBoundary(parent_mesh,
321 cylinder_surface_attributes);
322
323 char vishost[] = "localhost";
324 int visport = 19916;
325 socketstream cyl_sol_sock;
326 if (visualization)
327 {
328 cyl_sol_sock.open(vishost, visport);
329 cyl_sol_sock << "parallel " << num_procs << " " << myid << "\n";
330 cyl_sol_sock.precision(8);
331 cyl_sol_sock << "solution\n" << cylinder_submesh << temperature_cylinder_gf <<
332 "pause\n" << std::flush;
333 }
334 socketstream block_sol_sock;
335 if (visualization)
336 {
337 block_sol_sock.open(vishost, visport);
338 block_sol_sock << "parallel " << num_procs << " " << myid << "\n";
339 block_sol_sock.precision(8);
340 block_sol_sock << "solution\n" << block_submesh << temperature_block_gf <<
341 "pause\n" << std::flush;
342 }
343
344 // Create the transfer map needed in the time integration loop
345 auto temperature_block_to_cylinder_map = ParSubMesh::CreateTransferMap(
346 temperature_block_gf,
347 temperature_cylinder_gf);
348
349 real_t t = 0.0;
350 bool last_step = false;
351 for (int ti = 1; !last_step; ti++)
352 {
353 if (t + dt >= t_final - dt/2)
354 {
355 last_step = true;
356 }
357
358 // Advance the diffusion equation on the outer block to the next time step
359 d_ode_solver.Step(temperature_block, t, dt);
360 {
361 // Transfer the solution from the inner surface of the outer block to
362 // the cylinder outer surface to act as a boundary condition.
363 temperature_block_gf.SetFromTrueDofs(temperature_block);
364
365 temperature_block_to_cylinder_map.Transfer(temperature_block_gf,
366 temperature_cylinder_gf);
367
368 temperature_cylinder_gf.GetTrueDofs(temperature_cylinder);
369 }
370 // Advance the convection-diffusion equation on the outer block to the
371 // next time step
372 cd_ode_solver.Step(temperature_cylinder, t, dt);
373
374 if (last_step || (ti % vis_steps) == 0)
375 {
376 if (myid == 0)
377 {
378 out << "step " << ti << ", t = " << t << std::endl;
379 }
380
381 temperature_cylinder_gf.SetFromTrueDofs(temperature_cylinder);
382 temperature_block_gf.SetFromTrueDofs(temperature_block);
383
384 if (visualization)
385 {
386 cyl_sol_sock << "parallel " << num_procs << " " << myid << "\n";
387 cyl_sol_sock << "solution\n" << cylinder_submesh << temperature_cylinder_gf <<
388 std::flush;
389 block_sol_sock << "parallel " << num_procs << " " << myid << "\n";
390 block_sol_sock << "solution\n" << block_submesh << temperature_block_gf <<
391 std::flush;
392 }
393 }
394 }
395
396 return 0;
397}
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
Definition array.hpp:261
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T.
Definition array.hpp:269
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
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
Definition fespace.hpp:1313
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
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 Boundary Face Integrator. Assumes ownership of lfi.
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).
Pointer to an Operator of a specified type.
Definition handle.hpp:34
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition handle.hpp:145
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
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
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.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
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.
void SetFromTrueDofs(const Vector &tv) override
Set the GridFunction from the given true-dof vector.
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, const Array< int > &attr)
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
static ParTransferMap CreateTransferMap(const ParGridFunction &src, const ParGridFunction &dst)
Create a Transfer Map object.
static ParSubMesh CreateFromBoundary(const ParMesh &parent, Array< int > &boundary_attributes)
Create a surface ParSubMesh from it's parent.
Definition psubmesh.cpp:32
static ParSubMesh CreateFromDomain(const ParMesh &parent, Array< int > &domain_attributes)
Create a domain ParSubMesh from it's parent.
Definition psubmesh.cpp:26
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition ode.hpp:151
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition ode.cpp:68
void Step(Vector &x, real_t &t, real_t &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
Definition ode.cpp:76
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
Base abstract class for first order time dependent operators.
Definition operator.hpp:317
Base class for vector Coefficients that optionally depend on time and space.
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:247
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
const real_t alpha
Definition ex15.cpp:369
real_t kappa
Definition ex24.cpp:54
int main()
void velocity_profile(const Vector &c, Vector &q)
void velocity_profile(const Vector &c, Vector &q)
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
float real_t
Definition config.hpp:43
const char vishost[]
RefCoord t[3]