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;
83 ReducedSystemOperator *reduced_oper;
97 enum NonlinearSolverType
104 double visc,
double mu,
double K,
105 NonlinearSolverType nls_type);
111 virtual void ImplicitSolve(
const double dt,
const Vector &x,
Vector &k);
117 void InitSundialsJacSolver(SundialsJacSolver &sjsolv);
119 double ElasticEnergy(
const Vector &x)
const;
120 double KineticEnergy(
const Vector &v)
const;
123 virtual ~HyperelasticOperator();
130 class ReducedSystemOperator :
public Operator
144 void SetParameters(
double dt_,
const Vector *v_,
const Vector *x_);
152 virtual ~ReducedSystemOperator();
177 : M(), S(), H(), grad_H(), Jacobian(), J_solver() { }
183 M = &M_; S = &S_; H = &H_; J_solver = &solver;
193 int InitSystem(
void *sundials_mem);
194 int SetupSystem(
void *sundials_mem,
int conv_fail,
195 const Vector &y_pred,
const Vector &f_pred,
int &jac_cur,
197 int SolveSystem(
void *sundials_mem,
Vector &b,
const Vector &weight,
199 int FreeSystem(
void *sundials_mem);
205 class ElasticEnergyCoefficient :
public Coefficient
214 : model(m), x(x_) { }
216 virtual ~ElasticEnergyCoefficient() { }
225 bool init_vis =
false);
228 int main(
int argc,
char *argv[])
231 const char *mesh_file =
"../../data/beam-quad.mesh";
234 int ode_solver_type = 3;
235 double t_final = 300.0;
240 bool visualization =
true;
241 const char *nls =
"newton";
245 const double reltol = 1e-1, abstol = 1e-1;
248 args.
AddOption(&mesh_file,
"-m",
"--mesh",
249 "Mesh file to use.");
250 args.
AddOption(&ref_levels,
"-r",
"--refine",
251 "Number of times to refine the mesh uniformly.");
253 "Order (degree) of the finite elements.");
254 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
255 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
256 " 4 - CVODE implicit, approximate Jacobian,\n\t"
257 " 5 - CVODE implicit, specified Jacobian,\n\t"
258 " 6 - ARKODE implicit, approximate Jacobian,\n\t"
259 " 7 - ARKODE implicit, specified Jacobian,\n\t"
260 " 11 - Forward Euler, 12 - RK2,\n\t"
261 " 13 - RK3 SSP, 14 - RK4,\n\t"
262 " 15 - CVODE (adaptive order) explicit,\n\t"
263 " 16 - ARKODE default (4th order) explicit.");
264 args.
AddOption(&nls,
"-nls",
"--nonlinear-solver",
265 "Nonlinear systems solver: "
266 "\"newton\" (plain Newton) or \"kinsol\" (KINSOL).");
267 args.
AddOption(&t_final,
"-tf",
"--t-final",
268 "Final time; start time is 0.");
269 args.
AddOption(&dt,
"-dt",
"--time-step",
271 args.
AddOption(&visc,
"-v",
"--viscosity",
272 "Viscosity coefficient.");
273 args.
AddOption(&mu,
"-mu",
"--shear-modulus",
274 "Shear modulus in the Neo-Hookean hyperelastic model.");
275 args.
AddOption(&K,
"-K",
"--bulk-modulus",
276 "Bulk modulus in the Neo-Hookean hyperelastic model.");
277 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
278 "--no-visualization",
279 "Enable or disable GLVis visualization.");
280 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
281 "Visualize every n-th timestep.");
292 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
301 SundialsJacSolver *sjsolver = NULL;
302 switch (ode_solver_type)
313 if (ode_solver_type == 5)
315 sjsolver =
new SundialsJacSolver;
318 ode_solver = cvode;
break;
324 if (ode_solver_type == 7)
327 sjsolver =
new SundialsJacSolver;
330 ode_solver = arkode;
break;
333 case 12: ode_solver =
new RK2Solver(0.5);
break;
335 case 14: ode_solver =
new RK4Solver;
break;
340 ode_solver = cvode;
break;
345 ode_solver = arkode;
break;
351 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
356 map<string,HyperelasticOperator::NonlinearSolverType> nls_map;
357 nls_map[
"newton"] = HyperelasticOperator::NEWTON;
358 nls_map[
"kinsol"] = HyperelasticOperator::KINSOL;
359 if (nls_map.find(nls) == nls_map.end())
361 cout <<
"Unknown type of nonlinear solver: " << nls << endl;
368 for (
int lev = 0; lev < ref_levels; lev++)
383 cout <<
"Number of velocity/deformation unknowns: " << fe_size << endl;
386 fe_offset[1] = fe_size;
387 fe_offset[2] = 2*fe_size;
416 HyperelasticOperator oper(fespace, ess_bdr, visc, mu, K, nls_map[nls]);
421 char vishost[] =
"localhost";
423 vis_v.
open(vishost, visport);
426 visualize(vis_v, mesh, &x, &v,
"Velocity",
true);
427 vis_w.
open(vishost, visport);
430 oper.GetElasticEnergyDensity(x, w);
432 visualize(vis_w, mesh, &x, &w,
"Elastic energy density",
true);
438 cout <<
"initial elastic energy (EE) = " << ee0 << endl;
439 cout <<
"initial kinetic energy (KE) = " << ke0 << endl;
440 cout <<
"initial total energy (TE) = " << (ee0 + ke0) << endl;
444 ode_solver->
Init(oper);
448 bool last_step =
false;
449 for (
int ti = 1; !last_step; ti++)
451 double dt_real = min(dt, t_final - t);
453 ode_solver->
Step(vx, t, dt_real);
455 last_step = (t >= t_final - 1e-8*dt);
457 if (last_step || (ti % vis_steps) == 0)
462 cout <<
"step " << ti <<
", t = " << t <<
", EE = " << ee <<
", KE = "
463 << ke <<
", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
466 else if (arkode) { arkode->
PrintInfo(); }
474 oper.GetElasticEnergyDensity(x, w);
487 ofstream mesh_ofs(
"deformed.mesh");
488 mesh_ofs.precision(8);
489 mesh->
Print(mesh_ofs);
491 ofstream velo_ofs(
"velocity.sol");
492 velo_ofs.precision(8);
494 ofstream ee_ofs(
"elastic_energy.sol");
496 oper.GetElasticEnergyDensity(x, w);
510 GridFunction *field,
const char *field_name,
bool init_vis)
522 out <<
"solution\n" << *mesh << *field;
528 out <<
"window_size 800 800\n";
529 out <<
"window_title '" << field_name <<
"'\n";
536 out <<
"autoscale value\n";
543 ReducedSystemOperator::ReducedSystemOperator(
545 :
Operator(M_->Height()), M(M_), S(S_), H(H_), Jacobian(NULL),
546 dt(0.0), v(NULL), x(NULL), w(height), z(height)
549 void ReducedSystemOperator::SetParameters(
double dt_,
const Vector *v_,
552 dt = dt_; v = v_; x = x_;
565 Operator &ReducedSystemOperator::GetGradient(
const Vector &k)
const
568 Jacobian =
Add(1.0, M->SpMat(), dt, S->SpMat());
572 Jacobian->
Add(dt*dt, *grad_H);
576 ReducedSystemOperator::~ReducedSystemOperator()
582 int SundialsJacSolver::InitSystem(
void *sundials_mem)
585 HyperelasticOperator *he_oper;
588 he_oper =
dynamic_cast<HyperelasticOperator*
>(td_oper);
589 MFEM_VERIFY(he_oper,
"operator is not HyperelasticOperator");
594 he_oper->InitSundialsJacSolver(*
this);
598 int SundialsJacSolver::SetupSystem(
void *sundials_mem,
int conv_fail,
600 int &jac_cur,
Vector &v_temp1,
603 int sc = y_pred.
Size() / 2;
605 double dt = GetTimeStep(sundials_mem);
609 Jacobian =
Add(1.0, M->SpMat(), dt, S->SpMat());
610 grad_H =
dynamic_cast<SparseMatrix *
>(&H->GetGradient(x));
611 Jacobian->
Add(dt * dt, *grad_H);
613 J_solver->SetOperator(*Jacobian);
619 int SundialsJacSolver::SolveSystem(
void *sundials_mem,
Vector &b,
623 int sc = b.
Size() / 2;
628 double dt = GetTimeStep(sundials_mem);
631 grad_H->Mult(b_x, rhs);
633 M->AddMult(b_v, rhs);
635 J_solver->iterative_mode =
false;
636 J_solver->Mult(rhs, b_v);
643 int SundialsJacSolver::FreeSystem(
void *sundials_mem)
653 NonlinearSolverType nls_type)
655 M(&fespace), S(&fespace), H(&fespace),
656 viscosity(visc), z(height/2)
658 const double rel_tol = 1e-8;
659 const int skip_zero_entries = 0;
661 const double ref_density = 1.0;
664 M.Assemble(skip_zero_entries);
666 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
668 M.FormSystemMatrix(ess_tdof_list, tmp);
670 M_solver.iterative_mode =
false;
671 M_solver.SetRelTol(rel_tol);
672 M_solver.SetAbsTol(0.0);
673 M_solver.SetMaxIter(30);
674 M_solver.SetPrintLevel(0);
675 M_solver.SetPreconditioner(M_prec);
676 M_solver.SetOperator(M.SpMat());
680 H.SetEssentialTrueDofs(ess_tdof_list);
684 S.Assemble(skip_zero_entries);
685 S.FormSystemMatrix(ess_tdof_list, tmp);
687 reduced_oper =
new ReducedSystemOperator(&M, &S, &H);
689 #ifndef MFEM_USE_SUITESPARSE
703 if (nls_type == KINSOL)
707 newton_solver = kinsolver;
709 newton_solver->SetRelTol(rel_tol);
710 newton_solver->SetPrintLevel(0);
715 newton_solver->SetMaxIter(10);
716 newton_solver->SetRelTol(rel_tol);
717 newton_solver->SetPrintLevel(-1);
719 newton_solver->SetSolver(*J_solver);
720 newton_solver->iterative_mode =
false;
721 newton_solver->SetOperator(*reduced_oper);
734 if (viscosity != 0.0)
739 M_solver.Mult(z, dv_dt);
744 void HyperelasticOperator::ImplicitSolve(
const double dt,
759 reduced_oper->SetParameters(dt, &v, &x);
761 newton_solver->Mult(zero, dv_dt);
762 MFEM_VERIFY(newton_solver->GetConverged(),
763 "Nonlinear solver did not converge.");
765 cout <<
" num nonlin sol iters = " << newton_solver->GetNumIterations()
766 <<
", final norm = " << newton_solver->GetFinalNorm() <<
'\n';
768 add(v, dt, dv_dt, dx_dt);
771 void HyperelasticOperator::InitSundialsJacSolver(SundialsJacSolver &sjsolv)
773 sjsolv.SetOperators(M, S, H, *J_solver);
776 double HyperelasticOperator::ElasticEnergy(
const Vector &x)
const
778 return H.GetEnergy(x);
781 double HyperelasticOperator::KineticEnergy(
const Vector &v)
const
783 return 0.5*M.InnerProduct(v, v);
786 void HyperelasticOperator::GetElasticEnergyDensity(
789 ElasticEnergyCoefficient w_coeff(*model, x);
793 HyperelasticOperator::~HyperelasticOperator()
795 delete newton_solver;
806 model.SetTransformation(T);
809 return model.EvalW(J)/J.Det();
823 const double s = 0.1/64.;
826 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)
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)
virtual void Print(std::ostream &out=mfem::out) const
Conjugate gradient method.
Class for grid function - Vector with associated FE space.
Data type for scaled Jacobi-type smoother of sparse matrix.
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)
void SetLinearSolver(SundialsODELinearSolver &ls_spec)
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)
Direct sparse solver using UMFPACK.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
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.
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 int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
void PrintUsage(std::ostream &out) const
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
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)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
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
virtual void ProjectCoefficient(Coefficient &coeff)
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.
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.
void PrintInfo() const
Print ARKODE statistics.
Vector & GetBlock(int i)
Get the i-th vector in the block.
Arbitrary order "L2-conforming" discontinuous finite elements.
void Neg()
(*this) = -(*this)