MFEM  v4.6.0
Finite element discretization library
navier_cht.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // Overlapping Grids Miniapp: Conjugate heat transfer
14 // --------------------------------------------------
15 //
16 // This example code demonstrates use of MFEM to solve different physics in
17 // different domains using overlapping grids: A solid block with its base at a
18 // fixed temperature is cooled by incoming flow. The Fluid domain models the
19 // entire domain, minus the solid block, and the incompressible Navier-Stokes
20 // equations are solved on it:
21 //
22 // ________________________________________
23 // | |
24 // | FLUID DOMAIN |
25 // | |
26 // -->inflow | ______ | --> outflow
27 // (attr=1) | | | | (attr=2)
28 // |_______________| |_________________|
29 //
30 // Inhomogeneous Dirichlet conditions are imposed at inflow (attr=1) and
31 // homogeneous Dirichlet conditions are imposed on all surface (attr=3) except
32 // the outflow (attr=2) which has Neumann boundary conditions for velocity.
33 //
34 // In contrast to the Fluid domain, the Thermal domain includes the solid block,
35 // and the advection-diffusion equation is solved on it:
36 //
37 // dT/dt + u.grad T = kappa \nabla^2 T
38 //
39 // (attr=3)
40 // ________________________________________
41 // | |
42 // | THERMAL DOMAIN |
43 // (attr=1) | kappa1 |
44 // T=0 | ______ |
45 // | |kappa2| |
46 // |_______________|______|_________________|
47 // (attr=4) (attr=2) (attr=4)
48 // T=10
49 //
50 // Inhomogeneous boundary conditions (T=10) are imposed on the base of the solid
51 // block (attr=2) and homogeneous boundary conditions are imposed at the inflow
52 // region (attr=1). All other surfaces have Neumann condition.
53 //
54 // The one-sided coupling between the two domains is via transfer of the
55 // advection velocity (u) from fluid domain to thermal domain at each time step.
56 // mpirun -np 4 navier_cht -r1 3 -r2 2 -np1 2 -np2 2
57 
58 #include "mfem.hpp"
59 #include "navier_solver.hpp"
60 #include <fstream>
61 #include <iostream>
62 
63 using namespace std;
64 using namespace mfem;
65 using namespace navier;
66 
67 struct schwarz_common
68 {
69  // common
70  double dt = 2e-2;
71  double t_final = 250*dt;
72  // fluid
73  int fluid_order = 4;
74  double fluid_kin_vis = 0.001;
75  // solid
76  int solid_order = 4;
77  int ode_solver_type = 3;
78  double alpha = 1.0e-2;
79  double kappa = 0.5;
80 } schwarz;
81 
82 // Dirichlet conditions for velocity
83 void vel_dbc(const Vector &x, double t, Vector &u);
84 // solid conductivity
85 double kappa_fun(const Vector &x);
86 // initial condition for temperature
87 double temp_init(const Vector &x);
88 
89 class ConductionOperator : public TimeDependentOperator
90 {
91 protected:
92  ParFiniteElementSpace &fespace;
93  Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
94 
95  mutable ParBilinearForm *M;
96  ParBilinearForm *K;
97 
98  HypreParMatrix Mmat;
99  HypreParMatrix Kmat;
100  HypreParMatrix *T; // T = M + dt K
101  double current_dt;
102 
103  mutable CGSolver M_solver; // Krylov solver for inverting the mass matrix M
104  HypreSmoother M_prec; // Preconditioner for the mass matrix M
105 
106  CGSolver T_solver; // Implicit solver for T = M + dt K
107  HypreSmoother T_prec; // Preconditioner for the implicit solver
108 
109  double alpha, kappa, udir;
110 
111  mutable Vector z; // auxiliary vector
112 
113 public:
114  ConductionOperator(ParFiniteElementSpace &f, double alpha, double kappa,
116 
117  virtual void Mult(const Vector &u, Vector &du_dt) const;
118  /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
119  This is the only requirement for high-order SDIRK implicit integration.*/
120  virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
121 
122  /// Update the diffusion BilinearForm K using the given true-dof vector `u`.
123  void SetParameters(VectorGridFunctionCoefficient adv_gf_c);
124 
125  virtual ~ConductionOperator();
126 };
127 
128 void VisualizeField(socketstream &sock, const char *vishost, int visport,
129  ParGridFunction &gf, const char *title,
130  int x = 0, int y = 0, int w = 400, int h = 400,
131  bool vec = false);
132 
133 int main(int argc, char *argv[])
134 {
135  // Initialize MPI and HYPRE.
136  Mpi::Init(argc, argv);
137  int myid = Mpi::WorldRank();
138  Hypre::Init();
139 
140  // Parse command-line options.
141  int lim_meshes = 2; // should be greater than nmeshes
142  Array <const char *> mesh_file_list(lim_meshes);
143  Array <int> np_list(lim_meshes),
144  rs_levels(lim_meshes);
145  rs_levels = 0;
146  np_list = 1;
147  bool visualization = true;
148 
149  OptionsParser args(argc, argv);
150  args.AddOption(&np_list[0], "-np1", "--np1",
151  "number of MPI ranks for mesh 1");
152  args.AddOption(&np_list[1], "-np2", "--np2",
153  "number of MPI ranks for mesh 1");
154  args.AddOption(&rs_levels[0], "-r1", "--refine-serial 1",
155  "Number of times to refine the mesh 1 uniformly in serial.");
156  args.AddOption(&rs_levels[1], "-r2", "--refine-serial 2",
157  "Number of times to refine the mesh 2 uniformly in serial.");
158  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
159  "--no-visualization",
160  "Enable or disable GLVis visualization.");
161 
162  args.Parse();
163  if (!args.Good())
164  {
165  args.PrintUsage(cout);
166  return 1;
167  }
168  if (myid == 0)
169  {
170  args.PrintOptions(cout);
171  }
172 
173  const int nmeshes = 2;
174  mesh_file_list[0] = "fluid-cht.mesh";
175  mesh_file_list[1] = "solid-cht.mesh";
176 
177  // Setup MPI communicator for each mesh
178  MPI_Comm *comml = new MPI_Comm;
179  int color = 0;
180  int npsum = 0;
181  for (int i = 0; i < nmeshes; i++)
182  {
183  npsum += np_list[i];
184  if (myid < npsum) { color = i; break; }
185  }
186  MPI_Comm_split(MPI_COMM_WORLD, color, myid, comml);
187  int myidlocal, numproclocal;
188  MPI_Comm_rank(*comml, &myidlocal);
189  MPI_Comm_size(*comml, &numproclocal);
190 
191  Mesh *mesh = new Mesh(mesh_file_list[color], 1, 1);
192  int dim = mesh->Dimension();
193  mesh->SetCurvature(color == 0 ? schwarz.fluid_order : schwarz.solid_order);
194 
195  for (int lev = 0; lev < rs_levels[color]; lev++)
196  {
197  mesh->UniformRefinement();
198  }
199 
200 
201  if (color == 0 && myidlocal == 0)
202  {
203  std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
204  }
205 
206  // Setup ParMesh based on the communicator for each mesh
207  ParMesh *pmesh;
208  pmesh = new ParMesh(*comml, *mesh);
209  delete mesh;
210 
211  // Setup pointer for FESpaces, GridFunctions, and Solvers
212  H1_FECollection *fec_s = NULL; //FECollection for solid
213  ParFiniteElementSpace *fes_s = NULL; //FESpace for solid
214  ParFiniteElementSpace *adv_fes_s = NULL; //FESpace for advection in solid
215  ParGridFunction *u_gf = NULL; //Velocity solution on both meshes
216  ParGridFunction *t_gf = NULL; //Temperature solution
217  NavierSolver *flowsolver = NULL; //Fluid solver
218  ConductionOperator *coper = NULL; //Temperature solver
219  Vector t_tdof; //Temperature true-dof vector
220 
221  double t = 0,
222  dt = schwarz.dt,
223  t_final = schwarz.t_final;
224  bool last_step = false;
225 
226  // Setup flow solver on mesh for fluid
227  if (color == 0)
228  {
229  flowsolver = new NavierSolver(pmesh, schwarz.fluid_order,
230  schwarz.fluid_kin_vis);
231  flowsolver->EnablePA(true);
232  u_gf = flowsolver->GetCurrentVelocity();
233  Vector init_vel(dim);
234  init_vel = 0.;
235  VectorConstantCoefficient u_excoeff(init_vel);
236  u_gf->ProjectCoefficient(u_excoeff);
237 
238  // Dirichlet boundary conditions for fluid
239  Array<int> attr(pmesh->bdr_attributes.Max());
240  // Inlet is attribute 1.
241  attr[0] = 1;
242  // Walls is attribute 3.
243  attr[2] = 1;
244  flowsolver->AddVelDirichletBC(vel_dbc, attr);
245 
246  flowsolver->Setup(dt);
247  u_gf = flowsolver->GetCurrentVelocity();
248  }
249 
250  // Setup temperature solver for mesh on solid
251  ODESolver *ode_solver = NULL;
252  Vector vxyz;
253  if (color == 1)
254  {
255  switch (schwarz.ode_solver_type)
256  {
257  // Implicit L-stable methods
258  case 1: ode_solver = new BackwardEulerSolver; break;
259  case 2: ode_solver = new SDIRK23Solver(2); break;
260  case 3: ode_solver = new SDIRK33Solver; break;
261  // Explicit methods
262  case 11: ode_solver = new ForwardEulerSolver; break;
263  case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
264  case 13: ode_solver = new RK3SSPSolver; break;
265  case 14: ode_solver = new RK4Solver; break;
266  case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
267  // Implicit A-stable methods (not L-stable)
268  case 22: ode_solver = new ImplicitMidpointSolver; break;
269  case 23: ode_solver = new SDIRK23Solver; break;
270  case 24: ode_solver = new SDIRK34Solver; break;
271  default:
272  std::cout << "Unknown ODE solver type: " << schwarz.ode_solver_type << '\n';
273  delete mesh;
274  return 3;
275  }
276  fec_s = new H1_FECollection(schwarz.solid_order, dim);
277  fes_s = new ParFiniteElementSpace(pmesh, fec_s);
278  adv_fes_s = new ParFiniteElementSpace(pmesh, fec_s, 2);
279  t_gf = new ParGridFunction(fes_s);
280  u_gf = new ParGridFunction(adv_fes_s);
281 
283  t_gf->ProjectCoefficient(t_0);
284  t_gf->SetTrueVector();
285  t_gf->GetTrueDofs(t_tdof);
286 
287  // Create a list of points for the interior where the gridfunction will
288  // be interpolate from the fluid mesh
289  vxyz = *pmesh->GetNodes();
290  }
291 
292  // Setup FindPointsGSLIB. Note: we set it up with MPI_COMM_WORLD to enable
293  // communication between ParMesh for solid and fluid zones.
294  OversetFindPointsGSLIB finder(MPI_COMM_WORLD);
295  finder.Setup(*pmesh, color);
296 
297  // Tag each point to be found with the same id as the mesh
298  Array<unsigned int> color_array;
299  color_array.SetSize(vxyz.Size());
300  for (int i = 0; i < color_array.Size(); i++)
301  {
302  color_array[i] = (unsigned int)color;
303  }
304  Vector interp_vals(vxyz.Size());
305 
306  // Interpolate velocity solution on both meshes. Since the velocity solution
307  // does not exist on the temperature mesh, it just passes in a dummy
308  // gridfunction that is not used in any way on the fluid mesh.
309  finder.Interpolate(vxyz, color_array, *u_gf, interp_vals);
310 
311  // Transfer the interpolated solution to solid mesh and setup a coefficient.
313  if (color == 1)
314  {
315  *u_gf = interp_vals;
316  adv_gf_c.SetGridFunction(u_gf);
317  coper = new ConductionOperator(*fes_s, schwarz.alpha, schwarz.kappa,
318  adv_gf_c);
319  coper->SetParameters(adv_gf_c);
320  }
321 
322  // Visualize the solution.
323  char vishost[] = "localhost";
324  int visport = 19916;
325  socketstream vis_sol;
326  int Ww = 350, Wh = 350; // window size
327  int Wx = color*Ww+10, Wy = 0; // window position
328  if (visualization)
329  {
330  if (color == 0)
331  {
332  VisualizeField(vis_sol, vishost, visport, *u_gf,
333  "Velocity", Wx, Wy, Ww, Wh);
334  }
335  else
336  {
337  VisualizeField(vis_sol, vishost, visport, *t_gf,
338  "Temperature", Wx, Wy, Ww, Wh);
339  }
340  }
341 
342  if (ode_solver) { ode_solver->Init(*coper); }
343 
344  for (int step = 0; !last_step; ++step)
345  {
346  if (t + dt >= t_final - dt / 2)
347  {
348  last_step = true;
349  }
350 
351  double cfl;
352  if (flowsolver)
353  {
354  flowsolver->Step(t, dt, step);
355  cfl = flowsolver->ComputeCFL(*u_gf, dt);
356  }
357  if (ode_solver)
358  {
359  ode_solver->Step(t_tdof, t, dt);
360  t_gf->SetFromTrueDofs(t_tdof);
361  }
362  finder.Interpolate(vxyz, color_array, *u_gf, interp_vals);
363  if (color == 1)
364  {
365  *u_gf = interp_vals;
366  adv_gf_c.SetGridFunction(u_gf);
367  coper->SetParameters(adv_gf_c);
368  }
369 
370  if (visualization)
371  {
372  if (color == 0)
373  {
374  VisualizeField(vis_sol, vishost, visport, *u_gf,
375  "Velocity", Wx, Wy, Ww, Wh);
376  }
377  else
378  {
379  VisualizeField(vis_sol, vishost, visport, *t_gf,
380  "Temperature", Wx, Wy, Ww, Wh);
381  }
382  }
383 
384  if (color == 0 && myidlocal == 0)
385  {
386  printf("%11s %11s %11s\n", "Time", "dt", "CFL");
387  printf("%.5E %.5E %.5E\n", t, dt,cfl);
388  fflush(stdout);
389  }
390  }
391 
392  if (flowsolver) { flowsolver->PrintTimingData(); }
393 
394  finder.FreeData();
395  delete coper;
396  delete t_gf;
397  if (color == 1) { delete u_gf; }
398  delete adv_fes_s;
399  delete fes_s;
400  delete fec_s;
401  delete ode_solver;
402  delete flowsolver;
403  delete pmesh;
404  delete comml;
405 
406  return 0;
407 }
408 
409 ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
410  double kap,
412  : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
413  T(NULL), current_dt(0.0),
414  M_solver(f.GetComm()), T_solver(f.GetComm()), udir(10), z(height)
415 {
416  const double rel_tol = 1e-8;
417 
418  Array<int> ess_bdr(f.GetParMesh()->bdr_attributes.Max());
419  // Dirichlet boundary condition on inlet and isothermal section of wall.
420  ess_bdr = 0;
421  ess_bdr[0] = 1; // inlet
422  ess_bdr[1] = 1; // homogeneous isothermal section of bottom wall
423  ess_bdr[2] = 0; // top wall
424  ess_bdr[3] = 0; // inhomogeneous isothermal section of bottom wall
425  f.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
426 
427  M = new ParBilinearForm(&fespace);
428  M->AddDomainIntegrator(new MassIntegrator());
429  M->Assemble(0); // keep sparsity pattern of M and K the same
430  M->FormSystemMatrix(ess_tdof_list, Mmat);
431 
432  M_solver.iterative_mode = false;
433  M_solver.SetRelTol(rel_tol);
434  M_solver.SetAbsTol(0.0);
435  M_solver.SetMaxIter(100);
436  M_solver.SetPrintLevel(0);
437  M_prec.SetType(HypreSmoother::Jacobi);
438  M_solver.SetPreconditioner(M_prec);
439  M_solver.SetOperator(Mmat);
440 
441  alpha = al;
442  kappa = kap;
443 
444  T_solver.iterative_mode = false;
445  T_solver.SetRelTol(rel_tol);
446  T_solver.SetAbsTol(0.0);
447  T_solver.SetMaxIter(100);
448  T_solver.SetPrintLevel(0);
449  T_solver.SetPreconditioner(T_prec);
450 
451  SetParameters(adv_gf_c);
452 }
453 
454 void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
455 {
456  // Compute:
457  // du_dt = M^{-1}*-K(u)
458  // for du_dt
459 
460  Kmat.Mult(u, z);
461  z.Neg(); // z = -z
462  K->EliminateVDofsInRHS(ess_tdof_list, u, z);
463 
464  M_solver.Mult(z, du_dt);
465  du_dt.Print();
466  du_dt.SetSubVector(ess_tdof_list, 0.0);
467 }
468 
469 void ConductionOperator::ImplicitSolve(const double dt,
470  const Vector &u, Vector &du_dt)
471 {
472  // Solve the equation:
473  // du_dt = M^{-1}*[-K(u + dt*du_dt)]
474  // for du_dt
475  if (!T)
476  {
477  T = Add(1.0, Mmat, dt, Kmat);
478  current_dt = dt;
479  T_solver.SetOperator(*T);
480  }
481  MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
482  Kmat.Mult(u, z);
483  z.Neg();
484  K->EliminateVDofsInRHS(ess_tdof_list, u, z);
485 
486  T_solver.Mult(z, du_dt);
487  du_dt.SetSubVector(ess_tdof_list, 0.0);
488 }
489 
490 void ConductionOperator::SetParameters(VectorGridFunctionCoefficient adv_gf_c)
491 {
492  ParGridFunction u_alpha_gf(&fespace);
493  FunctionCoefficient kapfuncoef(kappa_fun);
494  u_alpha_gf.ProjectCoefficient(kapfuncoef);
495 
496  delete K;
497  K = new ParBilinearForm(&fespace);
498 
499  GridFunctionCoefficient u_coeff(&u_alpha_gf);
500 
501  K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
502  K->AddDomainIntegrator(new MixedDirectionalDerivativeIntegrator(adv_gf_c));
503  K->Assemble(0); // keep sparsity pattern of M and K the same
504  K->FormSystemMatrix(ess_tdof_list, Kmat);
505  delete T;
506  T = NULL; // re-compute T on the next ImplicitSolve
507 }
508 
509 ConductionOperator::~ConductionOperator()
510 {
511  delete T;
512  delete M;
513  delete K;
514 }
515 
516 void VisualizeField(socketstream &sock, const char *vishost, int visport,
517  ParGridFunction &gf, const char *title,
518  int x, int y, int w, int h, bool vec)
519 {
520  gf.HostRead();
521  ParMesh &pmesh = *gf.ParFESpace()->GetParMesh();
522  MPI_Comm comm = pmesh.GetComm();
523 
524  int num_procs, myid;
525  MPI_Comm_size(comm, &num_procs);
526  MPI_Comm_rank(comm, &myid);
527 
528  bool newly_opened = false;
529  int connection_failed;
530 
531  do
532  {
533  if (myid == 0)
534  {
535  if (!sock.is_open() || !sock)
536  {
537  sock.open(vishost, visport);
538  sock.precision(8);
539  newly_opened = true;
540  }
541  sock << "solution\n";
542  }
543 
544  pmesh.PrintAsOne(sock);
545  gf.SaveAsOne(sock);
546 
547  if (myid == 0 && newly_opened)
548  {
549  const char* keys = (gf.FESpace()->GetMesh()->Dimension() == 2)
550  ? "mAcRjlmm" : "mmaaAcl";
551 
552  sock << "window_title '" << title << "'\n"
553  << "window_geometry "
554  << x << " " << y << " " << w << " " << h << "\n"
555  << "keys " << keys;
556  if ( vec ) { sock << "vvv"; }
557  sock << std::endl;
558  }
559 
560  if (myid == 0)
561  {
562  connection_failed = !sock && !newly_opened;
563  }
564  MPI_Bcast(&connection_failed, 1, MPI_INT, 0, comm);
565  }
566  while (connection_failed);
567 }
568 
569 /// Fluid data
570 // Dirichlet conditions for velocity
571 void vel_dbc(const Vector &x, double t, Vector &u)
572 {
573  double xi = x(0);
574  double yi = x(1);
575 
576  u(0) = 0.;
577  u(1) = 0.;
578  if (std::fabs(xi+2.5)<1.e-5) { u(0) = 0.25*yi*(3-yi)/(1.5*1.5); }
579 }
580 
581 /// Solid data
582 // solid conductivity
583 double kappa_fun(const Vector &x)
584 {
585  return x(1) <= 1.0 && std::fabs(x(0)) < 0.5 ? 5.: 1.0;
586 }
587 
588 // initial temperature
589 double temp_init(const Vector &x)
590 {
591  double t_init = 1.0;
592  if (x(1) < 0.5)
593  {
594  t_init = 10*(std::exp(-x(1)*x(1)));
595  }
596  if (std::fabs(x(0)) >= 0.5)
597  {
598  double dx = std::fabs(x(0))-0.5;
599  t_init *= std::exp(-10*dx*dx);
600  }
601  return t_init;
602 }
void PrintTimingData()
Print timing summary of the solving routine.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
int visport
Conjugate gradient method.
Definition: solvers.hpp:493
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:935
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:457
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Base abstract class for first order time dependent operators.
Definition: operator.hpp:316
Vector coefficient that is constant in space and time.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:756
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
int Wx
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void Interpolate(const Vector &point_pos, Array< unsigned int > &point_id, const GridFunction &field_in, Vector &field_out, int point_pos_ordering=Ordering::byNODES)
Definition: gslib.cpp:1301
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:23
double kappa
Definition: ex24.cpp:54
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
STL namespace.
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:411
void Setup(double dt)
Initialize forms, solvers and preconditioners.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2320
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
double ComputeCFL(ParGridFunction &u, double dt)
Compute CFL.
OversetFindPointsGSLIB enables use of findpts for arbitrary number of overlapping grids...
Definition: gslib.hpp:231
virtual void FreeData()
Definition: gslib.cpp:272
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5635
void EnablePA(bool pa)
Enable partial assembly for every operator.
void AddVelDirichletBC(VectorCoefficient *coeff, Array< int > &attr)
Add a Dirichlet boundary condition to the velocity field.
Parallel smoothers in hypre.
Definition: hypre.hpp:971
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:163
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
bool is_open()
True if the socketstream is open, false otherwise.
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:150
void PrintAsOne(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4969
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:424
int Wh
int dim
Definition: ex24.cpp:53
void SetGridFunction(const GridFunction *gf)
Set the grid function for this coefficient. Also sets the Vector dimension to match that of the gf...
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:152
Vector coefficient defined by a vector GridFunction.
int Wy
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
int Ww
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
void Step(double &time, double dt, int cur_step, bool provisional=false)
Compute solution at the next time step t+dt.
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The classical forward Euler method.
Definition: ode.hpp:117
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
Transient incompressible Navier Stokes solver in a split scheme formulation.
void Setup(Mesh &m, const int meshid, GridFunction *gfmax=NULL, const double bb_t=0.1, const double newt_tol=1.0e-12, const int npt_max=256)
Definition: gslib.cpp:1135
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
double f(const Vector &p)