57 real_t r = sqrt(pow(x, 2.0) + pow(y, 2.0));
62 if (std::abs(r) >= 0.25 - 1e-8)
68 q(2) = A * exp(-(pow(x, 2.0) / 2.0 + pow(y, 2.0) / 2.0));
100 ess_tdofs_(ess_tdofs),
101 M_solver(fes.GetComm())
113 M.Reset(Mform.ParallelAssemble(),
true);
116 bform.AddBdrFaceIntegrator(
126 Kform.FormSystemMatrix(empty, K);
127 Mform.FormSystemMatrix(ess_tdofs_, M);
130 b = bform.ParallelAssemble();
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);
150 M_solver.Mult(t1, du_dt);
154 ~ConvectionDiffusionTDO()
204int main(
int argc,
char *argv[])
214 bool visualization =
true;
219 "Finite element order (polynomial degree).");
220 args.
AddOption(&t_final,
"-tf",
"--t-final",
221 "Final time; start time is 0.");
222 args.
AddOption(&dt,
"-dt",
"--time-step",
224 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
225 "--no-visualization",
226 "Enable or disable GLVis visualization.");
227 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
228 "Visualize every n-th timestep.");
231 Mesh *serial_mesh =
new Mesh(
"multidomain-hex.mesh");
243 cylinder_domain_attributes[0] = 1;
245 auto cylinder_submesh =
250 Array<int> inflow_attributes(cylinder_submesh.bdr_attributes.Max());
251 inflow_attributes = 0;
252 inflow_attributes[7] = 1;
255 cylinder_submesh.bdr_attributes.Max());
256 inner_cylinder_wall_attributes = 0;
257 inner_cylinder_wall_attributes[8] = 1;
262 Array<int> inflow_tdofs, interface_tdofs, ess_tdofs;
266 ess_tdofs.
Append(inflow_tdofs);
267 ess_tdofs.
Append(interface_tdofs);
270 ConvectionDiffusionTDO cd_tdo(fes_cylinder, ess_tdofs);
273 temperature_cylinder_gf = 0.0;
275 Vector temperature_cylinder;
276 temperature_cylinder_gf.
GetTrueDofs(temperature_cylinder);
279 cd_ode_solver.
Init(cd_tdo);
282 outer_domain_attributes[0] = 2;
285 outer_domain_attributes);
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;
297 block_submesh.bdr_attributes.Max());
298 outer_cylinder_wall_attributes = 0;
299 outer_cylinder_wall_attributes[8] = 1;
303 ConvectionDiffusionTDO d_tdo(fes_block, ess_tdofs, 0.0, 1.0);
306 temperature_block_gf = 0.0;
312 temperature_block_gf.
GetTrueDofs(temperature_block);
315 d_ode_solver.
Init(d_tdo);
318 cylinder_surface_attributes[0] = 9;
321 cylinder_surface_attributes);
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;
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;
346 temperature_block_gf,
347 temperature_cylinder_gf);
350 bool last_step =
false;
351 for (
int ti = 1; !last_step; ti++)
353 if (
t + dt >= t_final - dt/2)
359 d_ode_solver.
Step(temperature_block,
t, dt);
365 temperature_block_to_cylinder_map.Transfer(temperature_block_gf,
366 temperature_cylinder_gf);
368 temperature_cylinder_gf.
GetTrueDofs(temperature_cylinder);
372 cd_ode_solver.
Step(temperature_cylinder,
t, dt);
374 if (last_step || (ti % vis_steps) == 0)
378 out <<
"step " << ti <<
", t = " <<
t << std::endl;
386 cyl_sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
387 cyl_sol_sock <<
"solution\n" << cylinder_submesh << temperature_cylinder_gf <<
389 block_sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
390 block_sol_sock <<
"solution\n" << block_submesh << temperature_block_gf <<
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T.
Conjugate gradient method.
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)
Arbitrary order H1-conforming (continuous) finite elements.
Parallel smoothers in hypre.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
int Dimension() const
Dimension of the reference space used within the elements.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
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.
int height
Dimension of the output / number of rows in the matrix.
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...
Abstract parallel finite element space.
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
ParMesh * GetParMesh() const
Class for parallel grid function.
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 meshes.
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.
static ParSubMesh CreateFromDomain(const ParMesh &parent, Array< int > &domain_attributes)
Create a domain ParSubMesh from it's parent.
Third-order, strong stability preserving (SSP) Runge-Kutta method.
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
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].
Base abstract class for first order time dependent operators.
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x,...
Base class for vector Coefficients that optionally depend on time and space.
A general vector function coefficient.
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
void velocity_profile(const Vector &c, Vector &q)
void velocity_profile(const Vector &c, Vector &q)
real_t u(const Vector &xvec)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...