38#error This example requires that MFEM is built with MFEM_USE_PETSC=YES
91 virtual ~FE_Evolution() {
delete iJacobian;
delete rJacobian; }
117 MPI_Comm_size(pmesh->
GetComm(),&num_procs);
118 MPI_Comm_rank(pmesh->
GetComm(),&myid);
119 sout <<
"parallel " << num_procs <<
" " << myid <<
"\n";
120 sout <<
"solution\n" << *pmesh << *
u;
121 if (pause) { sout <<
"pause\n"; }
128 cout <<
"GLVis visualization paused."
129 <<
" Press space (in the GLVis window) to resume it.\n";
137int main(
int argc,
char *argv[])
147 const char *mesh_file =
"../../data/periodic-hexagon.mesh";
148 int ser_ref_levels = 2;
149 int par_ref_levels = 0;
154 const char *device_config =
"cpu";
155 int ode_solver_type = 4;
158 bool visualization =
true;
162 bool use_petsc =
true;
163 bool implicit =
false;
164 bool use_step =
true;
165 const char *petscrc_file =
"";
168 cout.precision(precision);
171 args.
AddOption(&mesh_file,
"-m",
"--mesh",
172 "Mesh file to use.");
174 "Problem setup to use. See options in velocity_function().");
175 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
176 "Number of times to refine the mesh uniformly in serial.");
177 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
178 "Number of times to refine the mesh uniformly in parallel.");
180 "Order (degree) of the finite elements.");
181 args.
AddOption(&pa,
"-pa",
"--partial-assembly",
"-no-pa",
182 "--no-partial-assembly",
"Enable Partial Assembly.");
183 args.
AddOption(&ea,
"-ea",
"--element-assembly",
"-no-ea",
184 "--no-element-assembly",
"Enable Element Assembly.");
185 args.
AddOption(&fa,
"-fa",
"--full-assembly",
"-no-fa",
186 "--no-full-assembly",
"Enable Full Assembly.");
187 args.
AddOption(&device_config,
"-d",
"--device",
188 "Device configuration string, see Device::Configure().");
189 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
190 "ODE solver: 1 - Forward Euler,\n\t"
191 " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
192 args.
AddOption(&t_final,
"-tf",
"--t-final",
193 "Final time; start time is 0.");
194 args.
AddOption(&dt,
"-dt",
"--time-step",
196 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
197 "--no-visualization",
198 "Enable or disable GLVis visualization.");
199 args.
AddOption(&visit,
"-visit",
"--visit-datafiles",
"-no-visit",
200 "--no-visit-datafiles",
201 "Save data files for VisIt (visit.llnl.gov) visualization.");
202 args.
AddOption(&binary,
"-binary",
"--binary-datafiles",
"-ascii",
204 "Use binary (Sidre) or ascii format for VisIt data files.");
205 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
206 "Visualize every n-th timestep.");
207 args.
AddOption(&use_petsc,
"-usepetsc",
"--usepetsc",
"-no-petsc",
209 "Use or not PETSc to solve the ODE system.");
210 args.
AddOption(&petscrc_file,
"-petscopts",
"--petscopts",
211 "PetscOptions file to use.");
212 args.
AddOption(&use_step,
"-usestep",
"--usestep",
"-no-step",
214 "Use the Step() or Run() method to solve the ODE system.");
215 args.
AddOption(&implicit,
"-implicit",
"--implicit",
"-no-implicit",
217 "Use or not an implicit method in PETSc to solve the ODE system.");
232 Device device(device_config);
233 if (myid == 0) { device.
Print(); }
237 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
244 UserMonitor *pmon = NULL;
247 switch (ode_solver_type)
250 case 2: ode_solver =
new RK2Solver(1.0);
break;
252 case 4: ode_solver =
new RK4Solver;
break;
253 case 6: ode_solver =
new RK6Solver;
break;
257 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
274 for (
int lev = 0; lev < ser_ref_levels; lev++)
289 for (
int lev = 0; lev < par_ref_levels; lev++)
302 cout <<
"Number of unknowns: " << global_vSize << endl;
338 b->AddBdrFaceIntegrator(
355 u->ProjectCoefficient(u0);
359 ostringstream mesh_name, sol_name;
360 mesh_name <<
"ex9-mesh." << setfill(
'0') << setw(6) << myid;
361 sol_name <<
"ex9-init." << setfill(
'0') << setw(6) << myid;
362 ofstream omesh(mesh_name.str().c_str());
363 omesh.precision(precision);
365 ofstream osol(sol_name.str().c_str());
366 osol.precision(precision);
380 MFEM_ABORT(
"Must build with MFEM_USE_SIDRE=YES for binary output.");
403 cout <<
"Unable to connect to GLVis server at "
405 visualization =
false;
408 cout <<
"GLVis visualization disabled.\n";
413 sout <<
"parallel " << num_procs <<
" " << myid <<
"\n";
414 sout.precision(precision);
415 sout <<
"solution\n" << *pmesh << *
u;
419 cout <<
"GLVis visualization paused."
420 <<
" Press space (in the GLVis window) to resume it.\n";
425 sout.precision(precision);
426 pmon =
new UserMonitor(sout,pmesh,
u,vis_steps);
432 FE_Evolution *adv =
new FE_Evolution(*m, *k, *B, implicit);
442 ode_solver->
Init(*adv);
450 for (
int ti = 0; !done; )
454 real_t dt_real = min(dt, t_final -
t);
455 ode_solver->
Step(*U,
t, dt_real);
458 done = (
t >= t_final - 1e-8*dt);
460 if (done || ti % vis_steps == 0)
464 cout <<
"time step: " << ti <<
", time: " <<
t << endl;
472 sout <<
"parallel " << num_procs <<
" " << myid <<
"\n";
473 sout <<
"solution\n" << *pmesh << *
u << flush;
485 else { ode_solver->
Run(*U,
t, dt, t_final); }
491 ostringstream sol_name;
492 sol_name <<
"ex9-final." << setfill(
'0') << setw(6) << myid;
493 ofstream osol(sol_name.str().c_str());
494 osol.precision(precision);
522 const Vector &b_,
bool M_in_lhs)
526 b(b_), comm(M_.ParFESpace()->GetComm()), M_solver(comm), z(height),
527 iJacobian(NULL), rJacobian(NULL)
542 M_solver.SetOperator(*M);
557 M_solver.SetPreconditioner(*M_prec);
558 M_solver.iterative_mode =
false;
559 M_solver.SetRelTol(1e-9);
560 M_solver.SetAbsTol(0.0);
561 M_solver.SetMaxIter(100);
562 M_solver.SetPrintLevel(0);
566void FE_Evolution::ExplicitMult(
const Vector &x,
Vector &y)
const
584void FE_Evolution::ImplicitMult(
const Vector &x,
const Vector &xp,
597void FE_Evolution::Mult(
const Vector &x,
Vector &y)
const
617 mfem_error(
"FE_Evolution::GetExplicitGradient(x): Capability not coded!");
636 mfem_error(
"FE_Evolution::GetImplicitGradient(x,xp,shift):"
637 " Capability not coded!");
649 for (
int i = 0; i <
dim; i++)
662 case 1: v(0) = 1.0;
break;
663 case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.);
break;
664 case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
676 case 1: v(0) = 1.0;
break;
677 case 2: v(0) = w*X(1); v(1) = -w*X(0);
break;
678 case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0;
break;
686 real_t d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
690 case 1: v(0) = 1.0;
break;
691 case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0);
break;
692 case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0;
break;
706 for (
int i = 0; i <
dim; i++)
720 return exp(-40.*pow(X(0)-0.5,2));
724 real_t rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
727 const real_t s = (1. + 0.25*cos(2*M_PI*X(2)));
731 return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
732 erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
738 real_t x_ = X(0), y_ = X(1), rho, phi;
741 return pow(sin(M_PI*rho),2)*sin(3*phi);
746 return sin(
f*X(0))*sin(
f*X(1));
@ GaussLobatto
Closed type.
Conjugate gradient method.
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
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)
virtual void Save()
Save the collection to disk.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
The classical forward Euler method.
A general function coefficient.
Wrapper for hypre's ParCSR matrix class.
Wrapper for hypre's parallel vector class.
Parallel smoothers in hypre.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Arbitrary order "L2-conforming" discontinuous finite elements.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
int Dimension() const
Dimension of the reference space used within the elements.
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.
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).
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
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].
virtual void Run(Vector &x, real_t &t, real_t &dt, real_t tf)
Perform time integration from time t [in] to time tf [in].
Pointer to an Operator of a specified type.
Jacobi smoothing for a given bilinear form (no matrix necessary).
Type
Enumeration defining IDs for some classes derived from Operator.
@ ANY_TYPE
ID for the base class Operator, i.e. any type.
@ PETSC_MATAIJ
ID for class PetscParMatrix, MATAIJ format.
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...
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
HYPRE_BigInt GlobalTrueVSize() const
Class for parallel grid function.
Class for parallel meshes.
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Abstract class for PETSc's ODE solvers.
virtual void Init(TimeDependentOperator &f_, enum PetscODESolver::Type type)
Initialize the ODE solver.
Wrapper for PETSc's matrix class.
Abstract class for monitoring PETSc's solvers.
virtual void MonitorSolution(PetscInt it, PetscReal norm, const Vector &x)
Monitor the solution vector x.
void SetMonitor(PetscSolverMonitor *ctx)
Sets user-defined monitoring routine.
Third-order, strong stability preserving (SSP) Runge-Kutta method.
The classical explicit forth-order Runge-Kutta method, RK4.
Data collection with Sidre routines following the Conduit mesh blueprint specification.
Base abstract class for first order time dependent operators.
virtual Operator & GetExplicitGradient(const Vector &x) const
Return an Operator representing dG/dx at the given point x and the currently set time.
virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const
Perform the action of the implicit part of the operator, F: y = F(x, k, t) where t is the current tim...
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,...
virtual void ExplicitMult(const Vector &x, Vector &y) const
Perform the action of the explicit part of the operator, G: y = G(x, t) where t is the current time.
virtual Operator & GetImplicitGradient(const Vector &x, const Vector &k, real_t shift) const
Return an Operator representing (dF/dk shift + dF/dx) at the given x, k, and the currently set time.
A general vector function coefficient.
int Size() const
Returns the size of the vector.
Data collection with VisIt I/O routines.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
real_t f(const Vector &p)
real_t u(const Vector &xvec)
void mfem_error(const char *msg)
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
void velocity_function(const Vector &x, Vector &v)
real_t inflow_function(const Vector &x)
real_t u0_function(const Vector &x)