37 #ifndef MFEM_USE_PETSC
38 #error This example requires that MFEM is built with MFEM_USE_PETSC=YES
83 virtual void ExplicitMult(
const Vector &x,
Vector &y)
const;
89 virtual ~FE_Evolution() {
delete iJacobian;
delete rJacobian; }
108 void MonitorSolution(PetscInt step, PetscReal norm,
const Vector &X)
115 MPI_Comm_size(pmesh->GetComm(),&num_procs);
116 MPI_Comm_rank(pmesh->GetComm(),&myid);
117 sout <<
"parallel " << num_procs <<
" " << myid <<
"\n";
118 sout <<
"solution\n" << *pmesh << *u;
119 if (pause) { sout <<
"pause\n"; }
126 cout <<
"GLVis visualization paused."
127 <<
" Press space (in the GLVis window) to resume it.\n";
135 int main(
int argc,
char *argv[])
139 MPI_Init(&argc, &argv);
140 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
141 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
145 const char *mesh_file =
"../../data/periodic-hexagon.mesh";
146 int ser_ref_levels = 2;
147 int par_ref_levels = 0;
149 int ode_solver_type = 4;
150 double t_final = 10.0;
152 bool visualization =
true;
156 bool use_petsc =
true;
157 bool implicit =
false;
158 bool use_step =
true;
159 const char *petscrc_file =
"";
162 cout.precision(precision);
165 args.
AddOption(&mesh_file,
"-m",
"--mesh",
166 "Mesh file to use.");
168 "Problem setup to use. See options in velocity_function().");
169 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
170 "Number of times to refine the mesh uniformly in serial.");
171 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
172 "Number of times to refine the mesh uniformly in parallel.");
174 "Order (degree) of the finite elements.");
175 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
176 "ODE solver: 1 - Forward Euler,\n\t"
177 " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
178 args.
AddOption(&t_final,
"-tf",
"--t-final",
179 "Final time; start time is 0.");
180 args.
AddOption(&dt,
"-dt",
"--time-step",
182 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
183 "--no-visualization",
184 "Enable or disable GLVis visualization.");
185 args.
AddOption(&visit,
"-visit",
"--visit-datafiles",
"-no-visit",
186 "--no-visit-datafiles",
187 "Save data files for VisIt (visit.llnl.gov) visualization.");
188 args.
AddOption(&binary,
"-binary",
"--binary-datafiles",
"-ascii",
190 "Use binary (Sidre) or ascii format for VisIt data files.");
191 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
192 "Visualize every n-th timestep.");
193 args.
AddOption(&use_petsc,
"-usepetsc",
"--usepetsc",
"-no-petsc",
195 "Use or not PETSc to solve the ODE system.");
196 args.
AddOption(&petscrc_file,
"-petscopts",
"--petscopts",
197 "PetscOptions file to use.");
198 args.
AddOption(&use_step,
"-usestep",
"--usestep",
"-no-step",
200 "Use the Step() or Run() method to solve the ODE system.");
201 args.
AddOption(&implicit,
"-implicit",
"--implicit",
"-no-implicit",
203 "Use or not an implicit method in PETSc to solve the ODE system.");
221 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
228 UserMonitor *pmon = NULL;
231 switch (ode_solver_type)
234 case 2: ode_solver =
new RK2Solver(1.0);
break;
236 case 4: ode_solver =
new RK4Solver;
break;
237 case 6: ode_solver =
new RK6Solver;
break;
241 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
251 PetscInitialize(NULL, NULL, petscrc_file, NULL);
259 for (
int lev = 0; lev < ser_ref_levels; lev++)
274 for (
int lev = 0; lev < par_ref_levels; lev++)
287 cout <<
"Number of unknowns: " << global_vSize << endl;
329 ostringstream mesh_name, sol_name;
330 mesh_name <<
"ex9-mesh." << setfill(
'0') << setw(6) << myid;
331 sol_name <<
"ex9-init." << setfill(
'0') << setw(6) << myid;
332 ofstream omesh(mesh_name.str().c_str());
333 omesh.precision(precision);
335 ofstream osol(sol_name.str().c_str());
336 osol.precision(precision);
347 #ifdef MFEM_USE_SIDRE
350 MFEM_ABORT(
"Must build with MFEM_USE_SIDRE=YES for binary output.");
367 char vishost[] =
"localhost";
369 sout.
open(vishost, visport);
373 cout <<
"Unable to connect to GLVis server at "
374 << vishost <<
':' << visport << endl;
375 visualization =
false;
378 cout <<
"GLVis visualization disabled.\n";
383 sout <<
"parallel " << num_procs <<
" " << myid <<
"\n";
384 sout.precision(precision);
385 sout <<
"solution\n" << *pmesh << *u;
389 cout <<
"GLVis visualization paused."
390 <<
" Press space (in the GLVis window) to resume it.\n";
395 sout.precision(precision);
396 pmon =
new UserMonitor(sout,pmesh,u,vis_steps);
402 FE_Evolution *adv =
new FE_Evolution(*M, *K, *B, implicit);
406 ode_solver->
Init(*adv);
413 for (
int ti = 0; !done; )
415 double dt_real = min(dt, t_final - t);
416 ode_solver->
Step(*U, t, dt_real);
419 done = (t >= t_final - 1e-8*dt);
421 if (done || ti % vis_steps == 0)
425 cout <<
"time step: " << ti <<
", time: " << t << endl;
434 sout <<
"parallel " << num_procs <<
" " << myid <<
"\n";
435 sout <<
"solution\n" << *pmesh << *u << flush;
447 else { ode_solver->
Run(*U, t, dt, t_final); }
453 ostringstream sol_name;
454 sol_name <<
"ex9-final." << setfill(
'0') << setw(6) << myid;
455 ofstream osol(sol_name.str().c_str());
456 osol.precision(precision);
478 if (use_petsc) { PetscFinalize(); }
487 const Vector &_b,
bool M_in_lhs)
491 M(_M), K(_K), b(_b), M_solver(M.GetComm()), z(_M.Height()),
492 iJacobian(NULL), rJacobian(NULL)
496 M_prec.SetType(HypreSmoother::Jacobi);
497 M_solver.SetPreconditioner(M_prec);
498 M_solver.SetOperator(M);
500 M_solver.iterative_mode =
false;
501 M_solver.SetRelTol(1e-9);
502 M_solver.SetAbsTol(0.0);
503 M_solver.SetMaxIter(100);
504 M_solver.SetPrintLevel(0);
509 void FE_Evolution::ExplicitMult(
const Vector &x,
Vector &y)
const
527 void FE_Evolution::ImplicitMult(
const Vector &x,
const Vector &xp,
549 Operator& FE_Evolution::GetExplicitGradient(
const Vector &x)
const
558 mfem_error(
"FE_Evolution::GetExplicitGradient(x): Capability not coded!");
575 mfem_error(
"FE_Evolution::GetImplicitGradient(x,xp,shift):"
576 " Capability not coded!");
588 for (
int i = 0; i <
dim; i++)
601 case 1: v(0) = 1.0;
break;
602 case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.);
break;
603 case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
612 const double w = M_PI/2;
615 case 1: v(0) = 1.0;
break;
616 case 2: v(0) = w*X(1); v(1) = -w*X(0);
break;
617 case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0;
break;
624 const double w = M_PI/2;
625 double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
629 case 1: v(0) = 1.0;
break;
630 case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0);
break;
631 case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0;
break;
645 for (
int i = 0; i <
dim; i++)
659 return exp(-40.*pow(X(0)-0.5,2));
663 double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
666 const double s = (1. + 0.25*cos(2*M_PI*X(2)));
670 return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
671 erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
677 double x_ = X(0), y_ = X(1), rho, phi;
680 return pow(sin(M_PI*rho),2)*sin(3*phi);
684 const double f = M_PI;
685 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.
Conjugate gradient method.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Wrapper for PETSc's matrix class.
Base abstract class for 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)
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]...
HYPRE_Int GlobalTrueVSize()
int Size() const
Returns the size of the vector.
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
virtual void Save(std::ostream &out) const
Abstract parallel finite element space.
void ProjectCoefficient(Coefficient &coeff)
void SetMonitor(PetscSolverMonitor *ctx)
Sets user-defined monitoring routine.
int main(int argc, char *argv[])
Abstract class for PETSc's ODE solvers.
Data collection with Sidre routines following the Conduit mesh blueprint specification.
Abstract class for monitoring PETSc's solvers.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Data collection with VisIt I/O routines.
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Parallel smoothers in hypre.
void PrintUsage(std::ostream &out) const
void SetTime(double t)
Set physical time (for time-dependent simulations)
Wrapper for hypre's parallel vector class.
The classical explicit forth-order Runge-Kutta method, RK4.
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...
void velocity_function(const Vector &x, Vector &v)
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)
void mfem_error(const char *msg)
Third-order, strong stability preserving (SSP) Runge-Kutta method.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
void PrintOptions(std::ostream &out) const
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
int open(const char hostname[], int port)
double u0_function(const Vector &x)
class for C-function coefficient
virtual void Run(Vector &x, double &t, double &dt, double tf)
Perform time integration from time t [in] to time tf [in].
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
Class for parallel grid function.
The classical forward Euler method.
Wrapper for hypre's ParCSR matrix class.
Class for parallel meshes.
double inflow_function(const Vector &x)
Arbitrary order "L2-conforming" discontinuous finite elements.
virtual void Print(std::ostream &out=std::cout) const