48 #ifndef MFEM_USE_SUNDIALS
49 #error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
55 class ReducedSystemOperator;
56 class SundialsJacSolver;
85 ReducedSystemOperator *reduced_oper;
99 enum NonlinearSolverType
106 double visc,
double mu,
double K,
107 NonlinearSolverType nls_type);
113 virtual void ImplicitSolve(
const double dt,
const Vector &x,
Vector &k);
119 void InitSundialsJacSolver(SundialsJacSolver &sjsolv);
126 virtual ~HyperelasticOperator();
133 class ReducedSystemOperator :
public Operator
149 void SetParameters(
double dt_,
const Vector *v_,
const Vector *x_);
157 virtual ~ReducedSystemOperator();
184 : M(), S(), H(), local_grad_H(), Jacobian(), J_solver() { }
191 M = &M_; S = &S_; H = &H_; J_solver = &solver;
192 ess_tdof_list = &ess_tdof_list_;
202 int InitSystem(
void *sundials_mem);
203 int SetupSystem(
void *sundials_mem,
int conv_fail,
204 const Vector &y_pred,
const Vector &f_pred,
int &jac_cur,
206 int SolveSystem(
void *sundials_mem,
Vector &b,
const Vector &weight,
208 int FreeSystem(
void *sundials_mem);
214 class ElasticEnergyCoefficient :
public Coefficient
223 : model(m), x(x_) { }
225 virtual ~ElasticEnergyCoefficient() { }
234 bool init_vis =
false);
237 int main(
int argc,
char *argv[])
241 MPI_Init(&argc, &argv);
242 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
243 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
246 const char *mesh_file =
"../../data/beam-quad.mesh";
247 int ser_ref_levels = 2;
248 int par_ref_levels = 0;
250 int ode_solver_type = 3;
251 double t_final = 300.0;
256 bool visualization =
true;
257 const char *nls =
"newton";
261 const double reltol = 1e-1, abstol = 1e-1;
264 args.
AddOption(&mesh_file,
"-m",
"--mesh",
265 "Mesh file to use.");
266 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
267 "Number of times to refine the mesh uniformly in serial.");
268 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
269 "Number of times to refine the mesh uniformly in parallel.");
271 "Order (degree) of the finite elements.");
272 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
273 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
274 " 4 - CVODE implicit, approximate Jacobian,\n\t"
275 " 5 - CVODE implicit, specified Jacobian,\n\t"
276 " 6 - ARKODE implicit, approximate Jacobian,\n\t"
277 " 7 - ARKODE implicit, specified Jacobian,\n\t"
278 " 11 - Forward Euler, 12 - RK2,\n\t"
279 " 13 - RK3 SSP, 14 - RK4,\n\t"
280 " 15 - CVODE (adaptive order) explicit,\n\t"
281 " 16 - ARKODE default (4th order) explicit.");
282 args.
AddOption(&nls,
"-nls",
"--nonlinear-solver",
283 "Nonlinear systems solver: "
284 "\"newton\" (plain Newton) or \"kinsol\" (KINSOL).");
285 args.
AddOption(&t_final,
"-tf",
"--t-final",
286 "Final time; start time is 0.");
287 args.
AddOption(&dt,
"-dt",
"--time-step",
289 args.
AddOption(&visc,
"-v",
"--viscosity",
290 "Viscosity coefficient.");
291 args.
AddOption(&mu,
"-mu",
"--shear-modulus",
292 "Shear modulus in the Neo-Hookean hyperelastic model.");
293 args.
AddOption(&K,
"-K",
"--bulk-modulus",
294 "Bulk modulus in the Neo-Hookean hyperelastic model.");
295 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
296 "--no-visualization",
297 "Enable or disable GLVis visualization.");
298 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
299 "Visualize every n-th timestep.");
318 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
327 SundialsJacSolver *sjsolver = NULL;
328 switch (ode_solver_type)
336 cvode =
new CVODESolver(MPI_COMM_WORLD, CV_BDF, CV_NEWTON);
339 if (ode_solver_type == 5)
341 sjsolver =
new SundialsJacSolver;
344 ode_solver = cvode;
break;
347 arkode =
new ARKODESolver(MPI_COMM_WORLD, ARKODESolver::IMPLICIT);
350 if (ode_solver_type == 7)
352 sjsolver =
new SundialsJacSolver;
355 ode_solver = arkode;
break;
358 case 12: ode_solver =
new RK2Solver(0.5);
break;
360 case 14: ode_solver =
new RK4Solver;
break;
362 cvode =
new CVODESolver(MPI_COMM_WORLD, CV_ADAMS, CV_FUNCTIONAL);
365 ode_solver = cvode;
break;
367 arkode =
new ARKODESolver(MPI_COMM_WORLD, ARKODESolver::EXPLICIT);
370 ode_solver = arkode;
break;
378 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
385 map<string,HyperelasticOperator::NonlinearSolverType> nls_map;
386 nls_map[
"newton"] = HyperelasticOperator::NEWTON;
387 nls_map[
"kinsol"] = HyperelasticOperator::KINSOL;
388 if (nls_map.find(nls) == nls_map.end())
392 cout <<
"Unknown type of nonlinear solver: " << nls << endl;
403 for (
int lev = 0; lev < ser_ref_levels; lev++)
413 for (
int lev = 0; lev < par_ref_levels; lev++)
430 cout <<
"Number of velocity/deformation unknowns: " << glob_size << endl;
435 true_offset[1] = true_size;
436 true_offset[2] = 2*true_size;
440 v_gf.
MakeTRef(&fespace, vx, true_offset[0]);
441 x_gf.
MakeTRef(&fespace, vx, true_offset[1]);
467 HyperelasticOperator oper(fespace, ess_bdr, visc, mu, K, nls_map[nls]);
472 char vishost[] =
"localhost";
474 vis_v.
open(vishost, visport);
476 visualize(vis_v, pmesh, &x_gf, &v_gf,
"Velocity",
true);
480 vis_w.
open(vishost, visport);
483 oper.GetElasticEnergyDensity(x_gf, w_gf);
485 visualize(vis_w, pmesh, &x_gf, &w_gf,
"Elastic energy density",
true);
489 double ee0 = oper.ElasticEnergy(x_gf);
490 double ke0 = oper.KineticEnergy(v_gf);
493 cout <<
"initial elastic energy (EE) = " << ee0 << endl;
494 cout <<
"initial kinetic energy (KE) = " << ke0 << endl;
495 cout <<
"initial total energy (TE) = " << (ee0 + ke0) << endl;
500 ode_solver->
Init(oper);
504 bool last_step =
false;
505 for (
int ti = 1; !last_step; ti++)
507 double dt_real = min(dt, t_final - t);
509 ode_solver->
Step(vx, t, dt_real);
511 last_step = (t >= t_final - 1e-8*dt);
513 if (last_step || (ti % vis_steps) == 0)
517 double ee = oper.ElasticEnergy(x_gf);
518 double ke = oper.KineticEnergy(v_gf);
522 cout <<
"step " << ti <<
", t = " << t <<
", EE = " << ee
523 <<
", KE = " << ke <<
", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
526 else if (arkode) { arkode->
PrintInfo(); }
534 oper.GetElasticEnergyDensity(x_gf, w_gf);
548 ostringstream mesh_name, velo_name, ee_name;
549 mesh_name <<
"deformed." << setfill(
'0') << setw(6) << myid;
550 velo_name <<
"velocity." << setfill(
'0') << setw(6) << myid;
551 ee_name <<
"elastic_energy." << setfill(
'0') << setw(6) << myid;
553 ofstream mesh_ofs(mesh_name.str().c_str());
554 mesh_ofs.precision(8);
555 pmesh->
Print(mesh_ofs);
557 ofstream velo_ofs(velo_name.str().c_str());
558 velo_ofs.precision(8);
560 ofstream ee_ofs(ee_name.str().c_str());
562 oper.GetElasticEnergyDensity(x_gf, w_gf);
590 out <<
"solution\n" << *mesh << *field;
596 out <<
"window_size 800 800\n";
597 out <<
"window_title '" << field_name <<
"'\n";
604 out <<
"autoscale value\n";
611 ReducedSystemOperator::ReducedSystemOperator(
614 :
Operator(M_->ParFESpace()->TrueVSize()), M(M_), S(S_), H(H_),
615 Jacobian(NULL), dt(0.0), v(NULL), x(NULL), w(height), z(height),
616 ess_tdof_list(ess_tdof_list_)
619 void ReducedSystemOperator::SetParameters(
double dt_,
const Vector *v_,
622 dt = dt_; v = v_; x = x_;
631 M->TrueAddMult(k, y);
632 S->TrueAddMult(w, y);
636 Operator &ReducedSystemOperator::GetGradient(
const Vector &k)
const
642 localJ->
Add(dt*dt, H->GetLocalGradient(z));
643 Jacobian = M->ParallelAssemble(localJ);
650 ReducedSystemOperator::~ReducedSystemOperator()
656 int SundialsJacSolver::InitSystem(
void *sundials_mem)
659 HyperelasticOperator *he_oper;
662 he_oper =
dynamic_cast<HyperelasticOperator*
>(td_oper);
663 MFEM_VERIFY(he_oper,
"operator is not HyperelasticOperator");
668 he_oper->InitSundialsJacSolver(*
this);
672 int SundialsJacSolver::SetupSystem(
void *sundials_mem,
int conv_fail,
674 int &jac_cur,
Vector &v_temp1,
677 int sc = y_pred.
Size() / 2;
679 double dt = GetTimeStep(sundials_mem);
684 local_grad_H = &H->GetLocalGradient(x);
685 localJ->
Add(dt*dt, *local_grad_H);
686 Jacobian = M->ParallelAssemble(localJ);
691 J_solver->SetOperator(*Jacobian);
697 int SundialsJacSolver::SolveSystem(
void *sundials_mem,
Vector &b,
701 int sc = b.
Size() / 2;
707 double dt = GetTimeStep(sundials_mem);
713 lb_x.Distribute(b_x);
714 local_grad_H->Mult(lb_x, lrhs);
715 lrhs.ParallelAssemble(rhs);
717 M->TrueAddMult(b_v, rhs);
718 rhs.SetSubVector(*ess_tdof_list, 0.0);
720 J_solver->iterative_mode =
false;
721 J_solver->Mult(rhs, b_v);
728 int SundialsJacSolver::FreeSystem(
void *sundials_mem)
738 NonlinearSolverType nls_type)
740 M(&fespace), S(&fespace), H(&fespace),
741 viscosity(visc), M_solver(f.GetComm()), z(height/2)
743 const double rel_tol = 1e-8;
744 const int skip_zero_entries = 0;
746 const double ref_density = 1.0;
749 M.Assemble(skip_zero_entries);
750 M.Finalize(skip_zero_entries);
751 Mmat = M.ParallelAssemble();
752 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
756 M_solver.iterative_mode =
false;
757 M_solver.SetRelTol(rel_tol);
758 M_solver.SetAbsTol(0.0);
759 M_solver.SetMaxIter(30);
760 M_solver.SetPrintLevel(0);
761 M_prec.SetType(HypreSmoother::Jacobi);
762 M_solver.SetPreconditioner(M_prec);
763 M_solver.SetOperator(*Mmat);
767 H.SetEssentialTrueDofs(ess_tdof_list);
771 S.Assemble(skip_zero_entries);
772 S.Finalize(skip_zero_entries);
774 reduced_oper =
new ReducedSystemOperator(&M, &S, &H, ess_tdof_list);
777 J_hypreSmoother->
SetType(HypreSmoother::l1Jacobi);
779 J_prec = J_hypreSmoother;
789 if (nls_type == KINSOL)
793 newton_solver = kinsolver;
795 newton_solver->SetRelTol(rel_tol);
796 newton_solver->SetPrintLevel(0);
802 newton_solver->SetRelTol(rel_tol);
803 newton_solver->SetPrintLevel(-1);
805 newton_solver->SetSolver(*J_solver);
806 newton_solver->iterative_mode =
false;
807 newton_solver->SetOperator(*reduced_oper);
820 if (viscosity != 0.0)
826 M_solver.Mult(z, dv_dt);
831 void HyperelasticOperator::ImplicitSolve(
const double dt,
846 reduced_oper->SetParameters(dt, &v, &x);
848 newton_solver->Mult(zero, dv_dt);
849 MFEM_VERIFY(newton_solver->GetConverged(),
850 "Nonlinear solver did not converge.");
852 if (fespace.GetMyRank() == 0)
854 cout <<
" num nonlin sol iters = " << newton_solver->GetNumIterations()
855 <<
", final norm = " << newton_solver->GetFinalNorm() <<
'\n';
858 add(v, dt, dv_dt, dx_dt);
861 void HyperelasticOperator::InitSundialsJacSolver(SundialsJacSolver &sjsolv)
863 sjsolv.SetOperators(M, S, H, *J_solver, ess_tdof_list);
866 double HyperelasticOperator::ElasticEnergy(
const ParGridFunction &x)
const
868 return H.GetEnergy(x);
871 double HyperelasticOperator::KineticEnergy(
const ParGridFunction &v)
const
873 double loc_energy = 0.5*M.InnerProduct(v, v);
875 MPI_Allreduce(&loc_energy, &energy, 1, MPI_DOUBLE, MPI_SUM,
880 void HyperelasticOperator::GetElasticEnergyDensity(
883 ElasticEnergyCoefficient w_coeff(*model, x);
887 HyperelasticOperator::~HyperelasticOperator()
889 delete newton_solver;
901 model.SetTransformation(T);
904 return model.EvalW(J)/J.Det();
918 const double s = 0.1/64.;
921 v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double t)
Evaluate the coefficient in the element described by T at the point ip at time t. ...
void Add(const int i, const int j, const double a)
void InitialDeformation(const Vector &x, Vector &y)
Conjugate gradient method.
Class for grid function - Vector with associated FE space.
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
A class to handle Vectors in a block fashion.
void SetMaxSetupCalls(int max_calls)
Set maximum number of nonlinear iterations without a Jacobian update.
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Subclass constant coefficient.
Base abstract class for time dependent operators.
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
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]...
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Data type dense matrix using column-major storage.
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
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Abstract parallel finite element space.
virtual void ProjectCoefficient(Coefficient &coeff)
int main(int argc, char *argv[])
Backward Euler ODE solver. L-stable.
double * GetData() const
Return a pointer to the beginning of the Vector data.
void add(const Vector &v1, const Vector &v2, Vector &v)
void SetSStolerances(double reltol, double abstol)
Specify the scalar relative and scalar absolute tolerances.
void InitialVelocity(const Vector &x, Vector &v)
void SetPositiveDiagonal(bool pos=true)
After computing l1-norms, replace them with their absolute values.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
void SetPrintLevel(int print_lvl)
Wrapper for SUNDIALS' KINSOL library – Nonlinear solvers.
Mesh * GetMesh() const
Returns the mesh.
Wrapper for SUNDIALS' CVODE library – Multi-step time integration.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
void SetMaxIter(int max_it)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
HYPRE_Int GlobalTrueVSize() const
void SetMaxStep(double dt_max)
Set the maximum time step of the linear multistep method.
Newton's method for solving F(x)=b for a given operator F.
virtual void Print(std::ostream &out=mfem::out) const
Parallel smoothers in hypre.
void PrintUsage(std::ostream &out) const
Wrapper for SUNDIALS' ARKODE library – Runge-Kutta time integration.
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
Set a custom Jacobian system solver for implicit methods.
int SpaceDimension() const
The classical explicit forth-order Runge-Kutta method, RK4.
void SetAbsTol(double atol)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void SetRelTol(double rtol)
Base class Coefficient that may optionally depend on time.
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)
Third-order, strong stability preserving (SSP) Runge-Kutta method.
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Abstract base class, wrapping the custom linear solvers interface in SUNDIALS' CVODE and ARKODE solve...
Implicit midpoint method. A-stable, not L-stable.
Class for integration point with weight.
void PrintOptions(std::ostream &out) const
Abstract class for hyperelastic models.
int open(const char hostname[], int port)
void GetNodes(Vector &node_coord) const
Arbitrary order H1-conforming (continuous) finite elements.
void PrintInfo() const
Print CVODE statistics.
Class for parallel grid function.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
The classical forward Euler method.
void SetMaxStep(double dt_max)
Set the maximum time step of the Runge-Kutta method.
Wrapper for hypre's ParCSR matrix class.
Class for parallel meshes.
void PrintInfo() const
Print ARKODE statistics.
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Arbitrary order "L2-conforming" discontinuous finite elements.