Commit 7c5abe96 authored by Lukas Riedel's avatar Lukas Riedel 🎧
Browse files

Merge branch 'refactoring/simulation'

parents a9aaf8fa 5d8c1ba9
......@@ -114,21 +114,43 @@ private:
using OWB = OutputWriterBase<Traits>;
Dune::ParameterTree& inifile;
typename Traits::GV& gv;
GFS& gfs;
Param& param;
U& u;
public:
/// Create the factory. References will be forwarded to the writer.
/** \param _inifile Parameter file parser
* \param _gv GridView
* \param _gfs GridFunctionSpace of the solution
* \param _param Parameter class
* \param _u Solution vector to write
*/
OutputWriterFactory (Dune::ParameterTree& _inifile, typename Traits::GV& _gv, GFS& _gfs, Param& _param, U& _u) :
inifile(_inifile),
gv(_gv),
gfs(_gfs),
param(_param),
u(_u)
{ }
/// Create a dummy OutputWriter
template<bool active = enabled>
static std::unique_ptr<OWB> create (Dune::ParameterTree& inifile, typename Traits::GV& gv, GFS& gfs, Param& param, U& unew, typename std::enable_if_t<!active,int> = 0)
std::enable_if_t<!active,std::unique_ptr<OWB>> create ()
{
return std::make_unique<OWB>();
}
/// Create a true OutputWriter
template<bool active = enabled>
static std::unique_ptr<OWB> create (Dune::ParameterTree& inifile, typename Traits::GV& gv, GFS& gfs, Param& param, U& unew, typename std::enable_if_t<active,int> = 0)
std::enable_if_t<active,std::unique_ptr<OWB>> create ()
{
std::unique_ptr<OWB> p;
p = std::make_unique<OutputWriter<Traits,GFS,Param,U>>(inifile,gv,gfs,param,unew);
p = std::make_unique<OutputWriter<Traits,GFS,Param,U>>(inifile,gv,gfs,param,u);
return p;
}
};
......
......@@ -44,7 +44,8 @@ Simulation<Traits>::Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid>
operator_setup();
// --- Utility Class Setup --- //
output = OutputWriterFactory::create(inifile,gv,*gfs,*param,*unew);
OutputWriterFactory output_fac(inifile,gv,*gfs,*param,*unew);
output = output_fac.create();
AdaptivityHandlerFactory adaptivity_fac(inifile,*grid);
adaptivity = adaptivity_fac.create();
......
......@@ -6,7 +6,6 @@
#include <dune/dorie/solver/util_grid_creator.hh>
#include <dune/dorie/solver/util_boundary_condition.hh>
#include <dune/dorie/solver/util_controller.hh>
#include <dune/dorie/solver/util_adaptivity.hh>
#include <dune/dorie/solver/util_h5tools.hh>
#include <dune/dorie/solver/richards_initial.hh>
#include <dune/dorie/solver/richards_boundary.hh>
......
#ifndef DUNE_DORIE_SOLVER_HH
#define DUNE_DORIE_SOLVER_HH
/// Assemble Operators and Routines for Solving. Then Solve the Problem
/** The actual program solving is executed after initialization of the necessary routines and operators.
* If the grid is adapted dynamically, the objects have to be re-built. Therefore, we use pointers for handling them.
* Calculation parameters, such as advancing time and catching errors, are handled by the Calculation Controller class.
* \tparam TRAITS Dorie::Traits type definitions
* \param grid Shared pointer to the grid
* \param gv GridView
* \param gfs Problem GridFunctionSpace
* \param afem Adaptivity FEM
* \param inifile Parameter File Parser
* \throw Dune::Exception Calculation does not run into fatal error (Errors in Newton Solver and Linear Solver can be handled)
*/
template<typename TRAITS, typename GridType, typename GV>
void RichardsSolver (GridType grid, GV& gv,
typename TRAITS::GFS& gfs,
typename TRAITS::AFEM& afem,
Dune::ParameterTree& inifile,
Dune::MPIHelper& helper)
{
Dune::Timer timer;
// Extract types for convenience
typedef typename TRAITS::RangeField RF;
enum { dim = TRAITS::dim };
// Problem space
typedef typename TRAITS::FEM FEM;
typedef typename TRAITS::CC CC;
typedef typename TRAITS::GFS GFS;
// Adaptivity space
typedef typename TRAITS::AFEM AFEM;
typedef typename TRAITS::AGFS AGFS;
// Read variables
const std::string filepath = inifile.get<std::string>("output.outputPath");
const std::string filename = inifile.get<std::string>("output.fileName");
const std::string gridType = inifile.get<std::string>("grid.gridType");
const int order = inifile.get<int>("grid.FEorder");
const bool useAdaptivity = inifile.get<bool>("adaptivity.useAdaptivity");
const bool asciiOutput = inifile.get<bool>("output.asciiVtk");
const int verbose = inifile.get<int>("output.verbose"); // we want verbosity per process
// Estimate matrix properties
RF entriesPerRow;
if(gridType=="rectangular"){
if(dim==2) entriesPerRow=5.0;
else if(dim==3) entriesPerRow=7.0;
}
else if(gridType=="gmsh"){
if(dim==2) entriesPerRow=13.0;
else if(dim==3) entriesPerRow=20.0;
}
gfs.name("matric_head");
CC cc;
// some output makes it look intelligent
gfs.update();
if(verbose>0){
if (gv.comm().rank()==0){
std::cout << "OPERATOR SETUP:" << std::endl;
std::cout << " Finite Element Order: " << order << std::endl;
std::cout << " Highest level on grid: " << grid->maxLevel() << std::endl;
}
auto DOFnumber = gfs.globalSize();
if(verbose>1)
std::cout << " Process " << gv.comm().rank() << ": number of DOF: " << DOFnumber << std::endl;
DOFnumber = gv.comm().sum(DOFnumber);
if (gv.comm().rank()==0)
std::cout << " Total number of DOF: " << DOFnumber << std::endl;
}
// --- Initialize correct parameterization class ---
typedef Dune::Dorie::InterpolatorBase<TRAITS,dim> INTERP;
typedef Dune::Dorie::InterpolatorFactory<TRAITS,dim> INTFAC;
INTFAC intfac(inifile,helper,verbose);
auto interp = intfac.create();
typedef Dune::Dorie::ParameterizationBase<TRAITS,INTERP> PARAM;
typedef Dune::Dorie::ParameterizationFactory<TRAITS,INTERP> PFAC;
PFAC pfac(inifile,helper,*interp,verbose);
auto param = pfac.create();
// --- Build Richards Classes ---
typedef Dune::Dorie::FlowBoundary<TRAITS> FBOUNDARY;
typedef Dune::Dorie::FlowSource<TRAITS,PARAM> FSOURCE;
FBOUNDARY fboundary(inifile);
FSOURCE fsource(inifile,*param);
// --- Set Up Grid Operators ---
typedef Dune::Dorie::Operator::RichardsDGSpatialOperator<TRAITS,PARAM,FBOUNDARY,FSOURCE,FEM,false> LOP;
typedef Dune::Dorie::Operator::RichardsDGTemporalOperator<TRAITS,PARAM,FEM,false> TLOP;
typedef Dune::PDELab::istl::BCRSMatrixBackend<RF> MBE;
typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,CC,CC> GO0;
typedef Dune::PDELab::GridOperator<GFS,GFS,TLOP,MBE,RF,RF,RF,CC,CC> GO1;
typedef Dune::PDELab::OneStepGridOperator<GO0,GO1> IGO;
typedef typename IGO::Traits::Domain U;
LOP lop(inifile,*param,gv,fboundary,fsource);
TLOP tlop(inifile,*param,gv);
MBE mbe_lop(entriesPerRow);
MBE mbe_tlop(1);
auto go0 = std::make_unique<GO0>(gfs,gfs,lop,mbe_lop);
auto go1 = std::make_unique<GO1>(gfs,gfs,tlop,mbe_tlop);
auto igo = std::make_unique<IGO>(*go0,*go1);
// --- Define Adaptivity Operator ---
typedef Dune::Dorie::FluxErrorEstimator<TRAITS,PARAM,AFEM> ESTLOP;
typedef Dune::PDELab::EmptyTransformation NoTrafo;
typedef typename Dune::PDELab::Backend::impl::BackendVectorSelector<AGFS,RF>::Type U0;
typedef Dune::PDELab::GridOperator<GFS,AGFS,ESTLOP,MBE,RF,RF,RF,NoTrafo,NoTrafo> ESTGO;
typedef Dune::Dorie::AdaptivityController<PARAM,AFEM,AGFS,ESTLOP,ESTGO,GV,GFS,GridType,U0,U> AC;
// --- Set Up Linear Solver ---
typedef Dune::PDELab::ISTLBackend_SEQ_SuperLU LS;
auto ls = std::make_unique<LS>(5000,0);
// --- Set Up Non-linear Newton Solver ---
typedef Dune::PDELab::Newton<IGO,LS,U> PDESOLVER;
auto pdesolver = std::make_unique<PDESOLVER>(*igo,*ls);
pdesolver->setParameters(inifile.sub("NewtonParameters"));
pdesolver->setVerbosityLevel(verbose);
// --- Set Up Time Step Scheme Operator ---
Dune::PDELab::Alexander2Parameter<RF> alex2;
typedef Dune::Dorie::OneStepMethod<RF,IGO,PDESOLVER,U,U> OSM;
auto osm = std::make_unique<OSM>(alex2,*igo,*pdesolver);
osm->setVerbosityLevel(verbose+1);
// --- Set Up Controller ---
Dune::Dorie::CalculationController<RF,FBOUNDARY> CCon(inifile,fboundary,helper);
// --- Build Solution vectors and interpolate with initial condition ---
U uold(gfs,0.0);
U unew(gfs,0.0);
Dune::Dorie::FlowInitial<TRAITS,PARAM> init(inifile,*param,gv);
Dune::PDELab::interpolate(init,gfs,uold);
unew = uold;
if(verbose>1 && gv.comm().rank()==0)
std::cout << "::: setup time " << std::setw(12) << timer.elapsed() << std::endl;
// --- Set Up VTK Plotting ---
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> UDGF;
UDGF udgf(gfs,unew);
auto udgf_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<UDGF>>(udgf,"head");
typedef Dune::Dorie::GradientFluxAdapter<GFS,U,PARAM> FluxDGF;
FluxDGF fluxdgf(gfs,unew,*param,udgf);
auto fluxdgf_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<FluxDGF>>(fluxdgf,"flux");
typedef Dune::Dorie::ConductivityAdapter<TRAITS,PARAM> ConDGF;
ConDGF condgf(gv,*param);
auto condgf_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ConDGF>>(condgf,"K_0");
typedef Dune::Dorie::WaterContentAdapter<TRAITS,PARAM,UDGF> WaterDGF;
WaterDGF waterdgf(gv,*param,udgf);
auto waterdgf_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<WaterDGF>>(waterdgf,"theta_w");
typedef Dune::Dorie::SaturationAdapter<TRAITS,PARAM,UDGF> SatDGF;
SatDGF satdgf(gv,*param,udgf);
auto satdgf_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<SatDGF>>(satdgf,"Theta");
Dune::VTKSequenceWriter<GV> vtkwriter(gv,filename,filepath,"vtk");
vtkwriter.addCellData(udgf_vtk);
vtkwriter.addCellData(fluxdgf_vtk);
vtkwriter.addCellData(condgf_vtk);
vtkwriter.addCellData(waterdgf_vtk);
vtkwriter.addCellData(satdgf_vtk);
try{ // VTK Plotter has very bad exception handling
RF t = CCon.getTime();
if(asciiOutput) vtkwriter.write(t,Dune::VTK::ascii);
else vtkwriter.write(t,Dune::VTK::base64);
}
catch(...){
DUNE_THROW(Dune::IOError,"Cannot write VTK output!");
}
// --- Calculate Solution ---
while(CCon.doStep())
{
timer.reset();
// get time values
RF t = CCon.getTime();
RF dt = CCon.getDT();
// limit newton solver
pdesolver->setMaxIterations(CCon.getIterations());
bool exception = false;
try
{
osm->apply(t,dt,uold,unew);
}
catch (Dune::PDELab::NewtonError &e){
exception = true;
std::cerr << " ERROR in Newton PDESolver: " << e << std::endl;
}
catch (Dune::ISTLError &e){
exception = true;
std::cerr << " ERROR in Linear Matrix Solver: " << e << std::endl;
}
catch (Dune::Exception &e){
std::cerr << "FATAL ERROR occured!" << std::endl;
DUNE_THROW(Dune::Exception,e);
}
// Controller receives information on success
if(CCon.validate(exception))
{
// print matrix statistics
if(verbose>2 && gv.comm().rank()==0){
typename IGO::Traits::Jacobian jac(*igo);
std::cout << jac.patternStatistics();
}
// --- Adaptive Grid Refinement ---
if (useAdaptivity)
{
AC ac(dim,inifile,*param,gv,afem);
ac.adapt(unew,gfs,grid);
Dune::Timer assemble_timer;
// Call load balance again to be sure
grid->loadBalance();
gfs.update();
if(verbose>1){
if (gv.comm().rank()==0){
std::cout << "OPERATOR SETUP:" << std::endl;
std::cout << " Highest level on grid: " << grid->maxLevel() << std::endl;
}
auto DOFnumber = gfs.globalSize();
if(verbose>1)
std::cout << " Process " << gv.comm().rank() << ": number of DOF: " << DOFnumber << std::endl;
DOFnumber = gv.comm().sum(DOFnumber);
if (gv.comm().rank()==0)
std::cout << " Total number of DOF: " << DOFnumber << std::endl;
}
// The grid has changed so we need to do the setup again
// Grid Operators
go0 = std::make_unique<GO0>(gfs,gfs,lop,mbe_lop);
go1 = std::make_unique<GO1>(gfs,gfs,tlop,mbe_tlop);
igo = std::make_unique<IGO>(*go0,*go1);
// Linear Solver
ls = std::make_unique<LS>(5000,0);
// Non-Linear Solver
pdesolver = std::make_unique<PDESOLVER>(*igo,*ls);
pdesolver->setParameters(inifile.sub("NewtonParameters"));
pdesolver->setVerbosityLevel(verbose);
// Time Step Scheme Operator
auto osm_result = osm->result();
osm = std::make_unique<OSM>(alex2,*igo,*pdesolver,osm_result);
osm->setStepNumber(osm_result.successful.timesteps+1);
osm->setVerbosityLevel(verbose+1);
if(verbose>1 && gv.comm().rank()==0)
std::cout << "::: setup time " << std::setw(12) << assemble_timer.elapsed() << std::endl;
} // adaptivity
// --- Write Solution ---
t += dt;
uold = unew;
if (gv.comm().rank()==0 && verbose>1){
std::cout << "::: entire time step " << std::setw(12) << timer.elapsed() << std::endl;
}
try{ // VTK Plotter has very bad exception handling
if(asciiOutput) vtkwriter.write(t,Dune::VTK::ascii);
else vtkwriter.write(t,Dune::VTK::base64);
}
catch(...){
DUNE_THROW(Dune::IOError,"Cannot write VTK output!");
}
}
}
}
#endif
#ifndef DUNE_DORIE_SOLVER_PREP_HH
#define DUNE_DORIE_SOLVER_PREP_HH
/// Define Grid Function Spaces for a simplex grid and call the Richards Solver Routine
/**
* \tparam RF Solution data type
* \tparam k Finite element order
* \param grid Shared pointer to the grid
* \param gv GridView
* \param inifile Parameter file reader
*/
template<typename RF, int k, typename GridType, typename GV>
void RichardsSolverSimplex (GridType grid, GV& gv, Dune::ParameterTree& inifile, Dune::MPIHelper& helper)
{
typedef typename GV::Grid::ctype DF;
enum { dim = GV::dimension };
// Create GFS for solving the problem
typedef Dune::PDELab::PkLocalFiniteElementMap<GV,DF,RF,k> FEM;
typedef Dune::PDELab::GridFunctionSpace<GV,FEM,
Dune::PDELab::NoConstraints,
Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed> >
GFS;
FEM fem(gv);
GFS gfs(gv,fem);
// Create GFS for adaptivity
typedef Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim> AFEM;
typedef Dune::PDELab::GridFunctionSpace<GV,AFEM,
Dune::PDELab::NoConstraints,
Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed> >
AGFS;
AFEM afem(Dune::GeometryType(Dune::GeometryType::simplex,dim));
typedef Dune::Dorie::Traits<GFS,AGFS> T;
RichardsSolver<T>(grid,gv,gfs,afem,inifile,helper);
}
/// Define Grid Function Spaces for a rectangular grid and call the Richards Solver Routine
/**
* \tparam RF Solution data type
* \tparam k Finite element order
* \param grid Shared pointer to the grid
* \param gv GridView
* \param inifile Parameter file reader
*/
template<typename RF, int k, typename GridType, typename GV>
void RichardsSolverRectangular (GridType grid, GV& gv, Dune::ParameterTree& inifile, Dune::MPIHelper& helper)
{
typedef typename GV::Grid::ctype DF;
enum { dim = GV::dimension };
// Create GFS for solving the problem
typedef Dune::PDELab::QkDGLocalFiniteElementMap<DF,RF,k,dim> FEM;
typedef Dune::PDELab::GridFunctionSpace<GV,FEM,
Dune::PDELab::NoConstraints,
Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed,
Dune::QkStuff::QkSize<k,dim>::value> >
GFS;
FEM fem;
GFS gfs(gv,fem);
// Create GFS for adaptivity
typedef Dune::PDELab::QkDGLocalFiniteElementMap<DF,RF,0,dim> AFEM;
typedef Dune::PDELab::GridFunctionSpace<GV,AFEM,
Dune::PDELab::NoConstraints,
Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed,
Dune::QkStuff::QkSize<0,dim>::value> >
AGFS;
AFEM afem;
typedef Dune::Dorie::Traits<GFS,AGFS> T;
RichardsSolver<T>(grid,gv,gfs,afem,inifile,helper);
}
#endif
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_UTIL_ADAPTIVITY_HH
#define DUNE_DORIE_UTIL_ADAPTIVITY_HH
namespace Dune {
namespace Dorie {
/// Class handling adaptive grid refinement with h-adaptivity
/** The jump in fluxes at a entity face is used as error estimator.
* The adaptation makes used of the PDELab-internal grid refinement functions.
*/
template<typename PARAM, typename P0FEM, typename P0GFS, typename ESTLOP, typename ESTGO, typename GV, typename GFS, typename Grid, typename U0, typename U>
class AdaptivityController
{
public:
const int maxLevel; //!< Maximum local grid refinement level
const int minLevel; //!< Minimum local grid refinement level
const std::string strategy; //!< Refinement strategy
const double alpha; //!< Refinement factor
const double beta; //!< Coarsening factor
const float adaptivityThreshold; //!< Global error threshold below which no refinement is applied
const int intorder; //!< Polynomial order of the local basis functions
int verbose; //!< Verbosity level
const int dim; //!< Spatial dimensions
PARAM& param; //!< Richards parameter class
const Dune::ParameterTree& inifile; //!< Parameter file
const GV& gv; //!< Grid View
const P0FEM& fem; //!< Finite Element Map for refinement space
/// Read adaptivity values from parameter file
/** \throw Dune::IOError Adaptive refinement strategy not supported
*/
AdaptivityController(const int dim_, const Dune::ParameterTree& inifile_, PARAM& param_, const GV& gv_, P0FEM& fem_)
: maxLevel(inifile_.get<int>("adaptivity.maxLevel")),
minLevel(inifile_.get<int>("adaptivity.minLevel")),
strategy(inifile_.get<std::string>("adaptivity.markingStrategy")),
alpha(inifile_.get<float>("adaptivity.refinementFraction")),
beta(inifile_.get<float>("adaptivity.coarseningFraction")),
adaptivityThreshold(inifile_.get<float>("adaptivity.threshold")),
intorder(inifile_.get<int>("grid.FEorder")),
verbose(inifile_.get<int>("output.verbose")),
dim(dim_), param(param_), inifile(inifile_), gv(gv_), fem(fem_)
{
if(strategy!="elementFraction" && strategy!="errorFraction" && strategy!="threshold" && strategy!="targetTolerance")
DUNE_THROW(Dune::IOError,"Valid values of adaptivity.markingStrategy: elementFraction, errorFraction, threshold, targetTolerance");
if(gv.comm().rank()!=0)
verbose=0;
}
/// Adapt the grid according to the estimated error and the strategy applied
/** \param u Solution vector
* \param gfs Problem Grid Function Space
* \param grid Shared pointer to the grid to be adapted
*/
void adapt(U& u, GFS& gfs, Grid grid)
{
Dune::Timer timer2;
float t2;
double maxeta(0.0);
int multiRefinementCounter = 0;
if(verbose>1){
std::cout << "GRID ADAPTION [" << strategy << "]:" << std::endl;
}
do { // Loop for multiple refinements
Dune::Timer timer3;
float t_setup,t_est,t_sqrt,t_strtgy,t_mark,t_adapt;
if(verbose>1 /*&& strategy == "targetTolerance"*/){
std::cout << " Refinement Step " << multiRefinementCounter << ": ";
}
double eta_alpha(0.0);
double eta_beta(0.0);
// Compute error estimate eta
P0GFS p0gfs(gv,fem);
ESTLOP estlop(inifile,param);
typename ESTGO::Traits::MatrixBackend mbe(0);
ESTGO estgo(gfs,p0gfs,estlop,mbe);
U0 eta(p0gfs,0.0);
t_setup = timer3.elapsed();
timer3.reset();
estgo.residual(u,eta);
t_est = timer3.elapsed();
timer3.reset();
using Dune::PDELab::Backend::native;
maxeta = native(eta)[0];
for (unsigned int i=0; i<eta.flatsize(); i++) {
native(eta)[i] = sqrt(native(eta)[i]); // eta contains squares
if (native(eta)[i] > maxeta) maxeta = native(eta)[i];
}
if (verbose>1) {
std::cout << "Max Local Error: " << maxeta << " " << std::endl;
if (maxeta < adaptivityThreshold) {
std::cout << " Max local error smaller than threshold " << adaptivityThreshold << ". Skipping grid refinement." << std::endl;
}
}
t_sqrt = timer3.elapsed();
timer3.reset();
// Stop adaptivity if target error is reached
if ((maxeta < alpha) && strategy == "targetTolerance") {
if (verbose>1)
std::cout << "Target error " << alpha << " reached!" << std::endl;
break;
}
// Apply marking strategy
if (strategy == "elementFraction")
Dune::PDELab::element_fraction( eta, alpha, beta, eta_alpha, eta_beta, verbose-1 );
else if (strategy == "errorFraction")
Dune::PDELab::error_fraction( eta, alpha, beta, eta_alpha, eta_beta, verbose-1 );
else{ //((strategy == "threshold") || (strategy == "targetTolerance")) {
eta_alpha = alpha;
eta_beta = beta; }
t_strtgy = timer3.elapsed();
timer3.reset();
// Skip refinement if threshold is met, but still allow coarsening
if (maxeta < adaptivityThreshold)
eta_alpha = maxeta + 1; // Only refine elements with error > maxerror + 1
Dune::PDELab::mark_grid( *grid, eta, eta_alpha, eta_beta, minLevel, maxLevel, verbose-1);
t_mark = timer3.elapsed();
timer3.reset();
Dune::PDELab::adapt_grid( *grid, gfs, u, intorder*2 );
t_adapt = timer3.elapsed();
timer3.reset();
multiRefinementCounter++;
if(verbose>3){
std::cout << "::: setup time" << std::setw(12) << t_setup << std::endl;
std::cout << "::: error estimation time" << std::setw(12) << t_est << std::endl;
std::cout << "::: unsquaring errors time" << std::setw(12) << t_sqrt << std::endl;
std::cout << "::: strategy application time" << std::setw(12) << t_strtgy << std::endl;
std::cout << "::: grid marking time" << std::setw(12) << t_mark << std::endl;
std::cout << "::: grid adaptation time" << std::setw(12) << t_adapt << std::endl;
}
} while (strategy == "targetTolerance");