69 class DG_Solver :
public Solver
81 BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
84 linear_solver.iterative_mode =
false;
85 linear_solver.SetRelTol(1e-9);
86 linear_solver.SetAbsTol(0.0);
87 linear_solver.SetMaxIter(100);
88 linear_solver.SetPrintLevel(0);
89 linear_solver.SetPreconditioner(prec);
92 void SetTimeStep(
double dt_)
103 linear_solver.SetOperator(A);
107 void SetOperator(
const Operator &op)
109 linear_solver.SetOperator(op);
114 linear_solver.Mult(x, y);
130 DG_Solver *dg_solver;
138 virtual void ImplicitSolve(
const double dt,
const Vector &x,
Vector &k);
144 int main(
int argc,
char *argv[])
148 const char *mesh_file =
"../data/periodic-hexagon.mesh";
154 const char *device_config =
"cpu";
155 int ode_solver_type = 4;
156 double t_final = 10.0;
158 bool visualization =
true;
160 bool paraview =
false;
165 cout.precision(precision);
168 args.
AddOption(&mesh_file,
"-m",
"--mesh",
169 "Mesh file to use.");
171 "Problem setup to use. See options in velocity_function().");
172 args.
AddOption(&ref_levels,
"-r",
"--refine",
173 "Number of times to refine the mesh uniformly.");
175 "Order (degree) of the finite elements.");
176 args.
AddOption(&pa,
"-pa",
"--partial-assembly",
"-no-pa",
177 "--no-partial-assembly",
"Enable Partial Assembly.");
178 args.
AddOption(&ea,
"-ea",
"--element-assembly",
"-no-ea",
179 "--no-element-assembly",
"Enable Element Assembly.");
180 args.
AddOption(&fa,
"-fa",
"--full-assembly",
"-no-fa",
181 "--no-full-assembly",
"Enable Full Assembly.");
182 args.
AddOption(&device_config,
"-d",
"--device",
183 "Device configuration string, see Device::Configure().");
184 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
185 "ODE solver: 1 - Forward Euler,\n\t"
186 " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6,\n\t"
187 " 11 - Backward Euler,\n\t"
188 " 12 - SDIRK23 (L-stable), 13 - SDIRK33,\n\t"
189 " 22 - Implicit Midpoint Method,\n\t"
190 " 23 - SDIRK23 (A-stable), 24 - SDIRK34");
191 args.
AddOption(&t_final,
"-tf",
"--t-final",
192 "Final time; start time is 0.");
193 args.
AddOption(&dt,
"-dt",
"--time-step",
195 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
196 "--no-visualization",
197 "Enable or disable GLVis visualization.");
198 args.
AddOption(&visit,
"-visit",
"--visit-datafiles",
"-no-visit",
199 "--no-visit-datafiles",
200 "Save data files for VisIt (visit.llnl.gov) visualization.");
201 args.
AddOption(¶view,
"-paraview",
"--paraview-datafiles",
"-no-paraview",
202 "--no-paraview-datafiles",
203 "Save data files for ParaView (paraview.org) visualization.");
204 args.
AddOption(&binary,
"-binary",
"--binary-datafiles",
"-ascii",
206 "Use binary (Sidre) or ascii format for VisIt data files.");
207 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
208 "Visualize every n-th timestep.");
217 Device device(device_config);
222 Mesh mesh(mesh_file, 1, 1);
228 switch (ode_solver_type)
232 case 2: ode_solver =
new RK2Solver(1.0);
break;
234 case 4: ode_solver =
new RK4Solver;
break;
235 case 6: ode_solver =
new RK6Solver;
break;
246 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
254 for (
int lev = 0; lev < ref_levels; lev++)
269 cout <<
"Number of unknowns: " << fes.
GetVSize() << endl;
296 constexpr
double alpha = -1.0;
321 ofstream omesh(
"ex9.mesh");
322 omesh.precision(precision);
324 ofstream osol(
"ex9-init.gf");
325 osol.precision(precision);
336 #ifdef MFEM_USE_SIDRE
339 MFEM_ABORT(
"Must build with MFEM_USE_SIDRE=YES for binary output.");
372 sout.
open(vishost, visport);
375 cout <<
"Unable to connect to GLVis server at "
376 << vishost <<
':' << visport << endl;
377 visualization =
false;
378 cout <<
"GLVis visualization disabled.\n";
382 sout.precision(precision);
383 sout <<
"solution\n" << mesh <<
u;
386 cout <<
"GLVis visualization paused."
387 <<
" Press space (in the GLVis window) to resume it.\n";
398 ode_solver->
Init(adv);
401 for (
int ti = 0; !done; )
403 double dt_real = min(dt, t_final - t);
404 ode_solver->
Step(u, t, dt_real);
407 done = (t >= t_final - 1e-8*dt);
409 if (done || ti % vis_steps == 0)
411 cout <<
"time step: " << ti <<
", time: " << t << endl;
415 sout <<
"solution\n" << mesh << u << flush;
437 ofstream osol(
"ex9-final.gf");
438 osol.precision(precision);
486 MFEM_VERIFY(dg_solver != NULL,
487 "Implicit time integration is not supported with partial assembly");
490 dg_solver->SetTimeStep(dt);
491 dg_solver->Mult(z, k);
507 for (
int i = 0; i <
dim; i++)
520 case 1: v(0) = 1.0;
break;
521 case 2: v(0) =
sqrt(2./3.); v(1) =
sqrt(1./3.);
break;
522 case 3: v(0) =
sqrt(3./6.); v(1) =
sqrt(2./6.); v(2) =
sqrt(1./6.);
531 const double w = M_PI/2;
534 case 1: v(0) = 1.0;
break;
535 case 2: v(0) = w*X(1); v(1) = -w*X(0);
break;
536 case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0;
break;
543 const double w = M_PI/2;
544 double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
548 case 1: v(0) = 1.0;
break;
549 case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0);
break;
550 case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0;
break;
564 for (
int i = 0; i <
dim; i++)
578 return exp(-40.*
pow(X(0)-0.5,2));
582 double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
585 const double s = (1. + 0.25*
cos(2*M_PI*X(2)));
589 return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
590 erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
596 double x_ = X(0), y_ = X(1), rho, phi;
603 const double f = M_PI;
604 return sin(f*X(0))*
sin(f*X(1));
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Conjugate gradient method.
Class for grid function - Vector with associated FE space.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Data type for scaled Jacobi-type smoother of sparse matrix.
void SetDataFormat(VTKFormat fmt)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Helper class for ParaView visualization data.
Base abstract class for first order time dependent operators.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
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 Size() const
Returns the size of the vector.
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
virtual void SetTime(const double t_)
Set the current time.
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
FDualNumber< tbase > exp(const FDualNumber< tbase > &f)
exp([dual number])
Backward Euler ODE solver. L-stable.
Data collection with Sidre routines following the Conduit mesh blueprint specification.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
virtual void Save()
Save the collection to disk.
Jacobi smoothing for a given bilinear form (no matrix necessary).
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
void SetHighOrderOutput(bool high_order_output_)
void PrintUsage(std::ostream &out) const
Print the usage message.
void SetTime(double t)
Set physical time (for time-dependent simulations)
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
A general vector function coefficient.
The classical explicit forth-order Runge-Kutta method, RK4.
void SetAbsTol(double atol)
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
void SetRelTol(double rtol)
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
void velocity_function(const Vector &x, Vector &v)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
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...
Third-order, strong stability preserving (SSP) Runge-Kutta method.
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...
Implicit midpoint method. A-stable, not L-stable.
FE_Evolution(FiniteElementSpace &vfes_, Operator &A_, SparseMatrix &Aflux_)
NURBSExtension * NURBSext
Optional NURBS mesh extension.
virtual void Print(std::ostream &os=mfem::out) const
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
void PrintOptions(std::ostream &out) const
Print the options.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
double u0_function(const Vector &x)
void SetLevelsOfDetail(int levels_of_detail_)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
A general function coefficient.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
virtual void Save() override
double u(const Vector &xvec)
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
The classical forward Euler method.
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
double inflow_function(const Vector &x)
Arbitrary order "L2-conforming" discontinuous finite elements.
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.