Merge branch 'development/dorie-v1.0' into 'master'

Various fixes to DG implementation and adaptation algorithms

Closes #26, #36, and #25

See merge request !14

dune/dorie/interface/adaptivity.hh

@@ -5,12 +5,13 @@ namespace Dune{
 namespace Dorie{
 
 /// Adaptivity base class. Does nothing.
-template<typename Traits, typename GFS, typename Param, typename U>
+template<typename Traits, typename GFS, typename Param, typename Boundary, typename U>
 class AdaptivityHandlerBase
 {
 private:
 	using Grid = typename Traits::Grid;
 	using GV = typename Traits::GV;
+	using RF = typename Traits::RF;
 
 public:
 	AdaptivityHandlerBase () = default;
@@ -19,12 +20,16 @@ public:
 	/// Do nothing.
 	/** \return Boolean if operators have to be updated
 	 */
-	virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) { return false; }
+	virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& uold, U& unew) { return false; }
 };
 
 /// Specialized SimulationBase class providing functions and members for adaptive grid refinement
-template<typename Traits, typename GFS, typename Param, typename U>
-class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,U>
+template<typename Traits,
+	typename GFS,
+	typename Param,
+	typename Boundary,
+	typename U>
+class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>
 {
 private:
 	static constexpr int order = Traits::fem_order;
@@ -37,7 +42,7 @@ private:
 	/// Adaptivity GFS
 	using AGFS = typename AGFSHelper::Type;
 	/// Error estimator local operator
-	using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,typename AGFSHelper::FEM>;
+	using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,Boundary>;
 	/// Empty constraints container
 	using NoTrafo = Dune::PDELab::EmptyTransformation;
 	/// Solution vector type
@@ -66,7 +71,7 @@ public:
 	 *  \param grid Grid to adapt (reference is not saved)
 	 */
 	AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
-		: AdaptivityHandlerBase<Traits,GFS,Param,U>(),
+		: AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>(),
 		inifile(_inifile),
 		maxLevel(inifile.get<int>("adaptivity.maxLevel")),
 		minLevel(inifile.get<int>("adaptivity.minLevel")),
@@ -95,7 +100,15 @@ public:
 	 *  \param unew Solution vector
 	 *  \return Boolean if grid operators have to be set up again (true if grid changed)
 	 */
-	bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) override
+	bool adapt_grid (
+		Grid& grid,
+		GV& gv,
+		GFS& gfs,
+		const Param& param,
+		const Boundary& boundary,
+		const RF time,
+		U& uold,
+		U& unew) override
 	{
 		Dune::Timer timer2;
 		float t2;
@@ -108,7 +121,7 @@ public:
 
 		do { // Loop for multiple refinements
 			Dune::Timer timer3;
-			float t_setup,t_est,t_sqrt,t_strtgy,t_mark,t_adapt;
+			float t_setup, t_sqrt, t_est, t_strtgy, t_mark, t_adapt;
 
 			if(verbose>1)
 				std::cout << "  Refinement Step " << multiRefinementCounter << ":  ";
@@ -116,9 +129,10 @@ public:
 			double eta_alpha(0.0);
 			double eta_beta(0.0);
 
-			// Compute error estimate eta
+			// set up helper GFS
 			AGFS p0gfs = AGFSHelper::create(gv);
-			ESTLOP estlop(inifile,param);
+			ESTLOP estlop(inifile, param, boundary);
+			estlop.setTime(time);
 			MBE mbe(0);
 			ESTGO estgo(gfs,p0gfs,estlop,mbe);
 			U0 eta(p0gfs,0.0);
@@ -126,17 +140,21 @@ public:
 			t_setup = timer3.elapsed();
 			timer3.reset();
 
+			// Compute error estimate eta
 			estgo.residual(uold,eta);
 
 			t_est = timer3.elapsed();
 			timer3.reset();
 
+			// unsquare errors
 			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];
-			}
+			auto& eta_nat = native(eta);
+			std::transform(eta_nat.begin(), eta_nat.end(), eta_nat.begin(),
+				[](const auto val) { return std::sqrt(val); }
+			);
+
+			// compute largest error
+			maxeta = eta.infinity_norm();
 			if (verbose>1)
 				std::cout << "Largest Local Error: " << maxeta << "   " << std::endl;
 
@@ -186,7 +204,7 @@ public:
 				std::cout << std::setprecision(4) << std::scientific;
 				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 << "::: unsquaring 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;
@@ -211,14 +229,14 @@ public:
  *  The constexpr boolean adapt_grid in Dune::Dorie::BaseTraits will
  *  define which create() function is enabled at compile time.
  */
-template<typename Traits, typename GFS, typename Param, typename U>
+template<typename Traits, typename GFS, typename Param, typename Boundary, typename U>
 struct AdaptivityHandlerFactory
 {
 private:
 	static constexpr bool enabled = Traits::adapt_grid;
 	using Grid = typename Traits::Grid;
 
-	using AHB = AdaptivityHandlerBase<Traits,GFS,Param,U>;
+	using AHB = AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>;
 
 	Dune::ParameterTree& inifile;
 	Grid& grid;
@@ -252,7 +270,7 @@ public:
 									std::is_same<Grid,Dune::UGGrid<3>>::value,
 									"Adaptivity is only implemented for UGGrid!");
 		std::unique_ptr<AHB> p;
-		p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,U>>(inifile,grid);
+		p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,Boundary,U>>(inifile,grid);
 		return p;
 	}
 };

dune/dorie/interface/simulation.cc

@@ -165,12 +165,15 @@ bool Simulation<Traits>::compute_time_step ()
 template<typename Traits>
 void Simulation<Traits>::run ()
 {
-	output->write_vtk_output(controller->getTime());
+	const auto t_start = controller->getTime();
+	output->write_vtk_output(t_start);
 	while(controller->doStep()) {
+		const auto time = controller->getTime();
+		const auto dt = controller->getDT();
 		if(!compute_time_step()){
 			continue;
 		}
-		if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *uold, *unew)){ // reset operators if grid changes
+		if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *fboundary, time+dt, *uold, *unew)){ // reset operators if grid changes
 			operator_setup();
 		}
 		output->write_vtk_output(controller->getTime());

dune/dorie/interface/simulation.hh

@@ -83,9 +83,9 @@ protected:
 	/// Factory for creating the output writer
 	using OutputWriterFactory = Dune::Dorie::OutputWriterFactory<Traits,GFS,Parameters,U>;
 	/// Grid Adaptivity Handler base class
-	using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<Traits,GFS,Parameters,U>;
+	using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<Traits,GFS,Parameters,FlowBoundary,U>;
 	/// Factory for creating the AdaptivityHandler
-	using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<Traits,GFS,Parameters,U>;
+	using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<Traits,GFS,Parameters,FlowBoundary,U>;
 
 	Dune::MPIHelper& helper;
 	std::shared_ptr<Grid> grid;

dune/dorie/solver/operator_DG.hh

@@ -328,17 +328,18 @@ public:
         const RF numFlux_n = satCond_n * ( - (gradu_n * normal) + penalty * jump);
         const RF numFlux   = omega_s * numFlux_s + omega_n * numFlux_n;
 
-        RF head_upwind;
-        if (numFlux > 0)
-          head_upwind = param.headMillerToRef(u_s, p_global_s);
-        else
-          head_upwind = param.headMillerToRef(u_n, p_global_n);
-
-        // compute saturation from matrix head
-        const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
-
-        // compute relative conductivity
-        const RF relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
+        // compute upwind conductivity factor
+        RF relCond_upwind;
+        if (numFlux > 0) {
+          const RF head_upwind = param.headMillerToRef(u_s, p_global_s);
+          const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
+          relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
+        }
+        else {
+          const RF head_upwind = param.headMillerToRef(u_n, p_global_n);
+          const RF saturation_upwind = param.saturation(head_upwind, p_global_n);
+          relCond_upwind = param.condFactor(saturation_upwind, p_global_n);
+        }
 
         // update residual (flux term)
         // diffusion term
@@ -473,10 +474,22 @@ public:
           // compute numerical flux estimation
           const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
 
+          RF head_upwind;
+          if (numFlux > 0)
+            head_upwind = param.headMillerToRef(u_s, p_global_s);
+          else
+            head_upwind = g; // Miller here?!
+
+          // compute saturation from matrix head
+          const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
+
+          // compute relative conductivity
+          const RF relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
+
           // switches between Neumann and Dirichlet bc
           // this choice is kept constant for one time step
           normal_flux = boundary.j(ig.intersection(),it.position(),time);
-          if(numFlux < normal_flux) {
+          if(relCond_upwind * numFlux < normal_flux) {
             bcType = BCType::dirichlet;
             bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
           }
@@ -556,7 +569,7 @@ public:
           // consistency term
           // + penalty term
           for (unsigned int i = 0; i<lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s,i, numFlux * phi_s[i] * factor);
+            r_s.accumulate(lfsv_s,i, relCond_upwind * numFlux * phi_s[i] * factor);
 
           // update residual (symmetry term)
           // (non-)symmetric IP term

dune/dorie/solver/operator_error_indicator.hh

@@ -13,24 +13,13 @@
 namespace Dune {
   namespace Dorie {
 
-    //! Weighting at grid interface for error estimation
-    struct EstimatorWeights
-    {
-      enum Type { weightsOn, //!< Use saturated conductivities at interface for weighting (default)
-                  weightsOff //!< No weighting
-                };
-    };
-
     /// Local operator for residual-based error estimation.
     /*
      * A call to residual() of a grid operator space will assemble
-     * the quantity \eta_T^2 for each cell. Note that the squares
-     * of the cell indicator \eta_T is stored. To compute the global
-     * error estimate sum up all values and take the square root.
+     * the quantity \eta_T for each cell.
      * Skeleton integral is mostly the same as in the DG operator.
-     *
      */
-    template<typename Traits, typename Parameter, typename FEM>
+    template<typename Traits, typename Parameter, typename Boundary>
     class FluxErrorEstimator
       : public Dune::PDELab::LocalOperatorDefaultFlags
     {
@@ -49,28 +38,27 @@ namespace Dune {
       typedef typename Traits::Domain             Domain;
       typedef typename Traits::IntersectionDomain ID;
 
-      typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
-      typedef Dune::PDELab::LocalBasisCache<LocalBasisType> Cache;
-
       // pattern assembly flags
       enum { doPatternVolume = false };
       enum { doPatternSkeleton = false };
 
       // residual assembly flags
       enum { doAlphaSkeleton  = true };
+      enum { doAlphaBoundary  = true };
 
-      Parameter& param;
-      EstimatorWeights::Type weights;
+      const Parameter& param;
+      const Boundary& boundary;
       const RF penalty_factor;
-      int intorderadd;
-      int quadrature_factor;
-
+      const int intorderadd;
+      const int quadrature_factor;
+      RF time;
 
-      FluxErrorEstimator (const Dune::ParameterTree& config, Parameter& param_,
-        EstimatorWeights::Type weights_ = EstimatorWeights::weightsOn,
+      FluxErrorEstimator (const Dune::ParameterTree& config,
+        const Parameter& param_, const Boundary& boundary_,
         int intorderadd_ = 2, int quadrature_factor_ = 2)
-      : param(param_), penalty_factor(config.get<RF>("dg.penaltyFactor")),
-      intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
+      : param(param_), boundary(boundary_),
+        penalty_factor(config.get<RF>("dg.penaltyFactor")),
+        intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
       { }
 
       /// skeleton integral depending on test and ansatz functions
@@ -102,8 +90,7 @@ namespace Dune {
         RF sum(0.0);
 
         // loop over quadrature points and integrate normal flux
-        const auto quadrature_rule = quadratureRule(gtface,intorder);
-        for (const auto& it : quadrature_rule)
+        for (const auto& it : quadratureRule(gtface,intorder))
         {
           // retrieve unit normal vector
           const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it.position());
@@ -175,14 +162,6 @@ namespace Dune {
           // compute jump in solution
           const RF h_jump = u_s - u_n;
 
-          // compute weights
-          RF omega_s = 0.5, omega_n = 0.5;
-          if(weights == EstimatorWeights::weightsOn)
-          {
-            omega_s = satCond_n / (satCond_s + satCond_n);
-            omega_n = satCond_s / (satCond_s + satCond_n);
-          }
-
           // apply gravity vector
           gradu_s -= param.gravity();
           gradu_n -= param.gravity();
@@ -190,34 +169,139 @@ namespace Dune {
           // penalty term
           const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
 
-          // compute numerical flux
-          const RF numFlux_s = satCond_s * ( - (gradu_s * n_F) + penalty * h_jump);
-          const RF numFlux_n = satCond_n * ( - (gradu_n * n_F) + penalty * h_jump);
-
-          // compute jump in physical flux
-          const RF jump = omega_s * condFactor_s * numFlux_s - omega_n * condFactor_n * numFlux_n;
-
           // integration factor
           const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
 
-          // add square of error
-          sum += jump*jump*factor*factor;
+          const RF grad_jump = (satCond_s * condFactor_s * (n_F * gradu_s))
+            - (satCond_n * condFactor_n * (n_F * gradu_n));
+
+          const RF gamma = 2.0 * satCond_s * condFactor_s * satCond_n * condFactor_n
+            / ( satCond_s * condFactor_s + satCond_n * condFactor_n );
+          const RF head_jump = gamma * penalty * h_jump;
+
+          const RF total_jump = grad_jump/2.0 + head_jump;
+
+          sum += total_jump * total_jump * factor;
         }
 
-        // correct for the number of neighbors
-        sum *= 1. / quadrature_rule.size();
+        DF h_T_s = diameter(ig.inside().geometry());
+        DF h_T_n = diameter(ig.outside().geometry());
 
-        // what is this variable for?
-        DF h_T = 1;//std::max( diameter(ig.inside().geometry()),
-                  //         diameter(ig.outside().geometry()) );
+        DF C_P_T = 1.0 / M_PI;
+        DF C_F_T_s = h_T_s * ig.geometry().volume() / ig.inside().geometry().volume();
+        C_F_T_s *= C_P_T * (2.0 / dim + C_P_T);
+        DF C_F_T_n = h_T_n * ig.geometry().volume() / ig.outside().geometry().volume();
+        C_F_T_n *= C_P_T * (2.0 / dim + C_P_T);
 
         // accumulate indicator
-        r_s.accumulate(lfsv_s,0,h_T*sum);
-        r_n.accumulate(lfsv_n,0,h_T*sum);
+        r_s.accumulate(lfsv_s,0, h_T_s * C_F_T_s * sum );
+        r_n.accumulate(lfsv_n,0, h_T_n * C_F_T_n * sum );
+      }
+
+    template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
+    void alpha_boundary (const IG& ig,
+        const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
+        R& r_s) const
+    {
+      // get polynomial degree
+      const int degree = lfsu_s.finiteElement().localBasis().order();
+      const int intorder = intorderadd + quadrature_factor * degree;
+
+      // geometric factor of penalty
+      const RF h_F = ig.inside().geometry().volume() / ig.geometry().volume(); // Houston!
+
+      // get element geometry
+      auto gtface = ig.geometry();
+
+      // container for sum of errors
+      RF sum(0.0);
+
+      // loop over quadrature points and integrate normal flux
+      for (const auto& it : quadratureRule(gtface,intorder))
+      {
+        // check boundary condition type
+        const auto bc = boundary.bc(ig.intersection(), it.position(), time);
+        if (!BoundaryCondition::isDirichlet(bc)){
+          continue;
+        }
+
+        // position of quadrature point in local coordinates of elements
+        const Domain local_s = ig.geometryInInside().global(it.position());
+
+        // retrieve unit normal vector
+        const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it.position());
+
+        // evaluate basis functions
+        std::vector<Scalar> phi_s(lfsu_s.size());
+        lfsu_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
+
+        // evaluate u
+        RF u_s = 0.;
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          u_s += x_s(lfsu_s,i) * phi_s[i];
+
+        // evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
+        typedef typename LFSU::Traits::FiniteElementType::Traits
+          ::LocalBasisType::Traits::JacobianType JacobianType;
+        std::vector<JacobianType> gradphi_s(lfsu_s.size());
+        lfsu_s.finiteElement().localBasis().evaluateJacobian(local_s,gradphi_s);
+
+        // transform gradients to real element
+        Dune::FieldMatrix<DF,dim,dim> jac;
+        std::vector<Vector> tgradphi_s(lfsu_s.size());
+        jac = ig.inside().geometry().jacobianInverseTransposed(local_s);
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          jac.mv(gradphi_s[i][0],tgradphi_s[i]);
+
+        // compute gradient of u
+        Vector gradu_s(0.);
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
+
+        // offset for parameter queries
+        const RF eps = 1e-9;
+
+        // parameter queries: inside
+        Domain xGlobal_s = ig.inside().geometry().global(local_s);
+        xGlobal_s = xGlobal_s.axpy(-eps,n_F);
+
+        const RF head_s = param.headMillerToRef(u_s,xGlobal_s);
+        const RF satCond_s = param.condRefToMiller(param.K(xGlobal_s),xGlobal_s);
+        const RF saturation_s = param.saturation(head_s,xGlobal_s);
+        const RF condFactor_s = param.condFactor(saturation_s,xGlobal_s);
+
+        // boundary condition at position
+        const RF g = boundary.g(ig.intersection(), it.position(), time);
+
+        // compute jump in solution
+        const RF h_jump = u_s - g;
+
+        // apply gravity vector
+        gradu_s -= param.gravity();
+
+        // penalty term
+        const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
+
+        // integration factor
+        const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
+
+        const RF head_jump = satCond_s * condFactor_s * penalty * h_jump;
+
+        sum += head_jump * head_jump * factor;
       }
 
+      DF h_T_s = diameter(ig.inside().geometry());
+
+      DF C_P_T = 1.0 / M_PI;
+      DF C_F_T_s = h_T_s * ig.geometry().volume() / ig.inside().geometry().volume();
+      C_F_T_s *= C_P_T * (2.0 / dim + C_P_T);
+
+      // accumulate indicator
+      r_s.accumulate(lfsv_s,0, h_T_s * C_F_T_s * sum );
+    }
+
     private:
-      // probably redundant
+      //! compute diameter of a grid cell
       template<class GEO>
       typename GEO::ctype diameter (const GEO& geo) const
       {
@@ -237,6 +321,13 @@ namespace Dune {
         return hmax;
       }
 
+    public:
+      //! set operator time
+      void setTime (RF t)
+      {
+        time = t;
+      }
+
     };
 
   }

testing/CMakeLists.txt

@@ -11,18 +11,12 @@ endfunction()
 
 # dorie run
 dorie_add_system_test(dorie ode_homogeneous_1.mini)
-dorie_add_system_test(dorie ode_homogeneous_1_adapt.mini)
 dorie_add_system_test(dorie ode_homogeneous_2.mini)
-dorie_add_system_test(dorie ode_homogeneous_2_adapt.mini)
 dorie_add_system_test(dorie ode_homogeneous_3.mini)
-dorie_add_system_test(dorie ode_homogeneous_3_adapt.mini)
 
 dorie_add_system_test(dorie ode_heterogeneous_1.mini)
-dorie_add_system_test(dorie ode_heterogeneous_1_adapt.mini)
 dorie_add_system_test(dorie ode_heterogeneous_2.mini)
-dorie_add_system_test(dorie ode_heterogeneous_2_adapt.mini)
 dorie_add_system_test(dorie ode_heterogeneous_3.mini)
-dorie_add_system_test(dorie ode_heterogeneous_3_adapt.mini)
 
 dorie_add_system_test(dorie ode_evaporation.mini)
 dorie_add_system_test(dorie ode_linear_interpolator.mini)

testing/ode_heterogeneous_1_adapt.mini.in

-include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
-
-__name = ode_heterogeneous_1_adapt
-_test_command = run
-_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
-_evaluation = ode
-
-output.fileName = ode_heterogeneous_1_adapt | unique
-output.outputPath = ode_heterogeneous_1_adapt | unique
-output.verbose = 0
-
-boundary.file = "{_asset_path}/bcs/infiltration_2d.dat"
-
-parameters.arrayFile = "{_asset_path}/parfields/layered.h5"
-
-time.end = 1E7
-time.maxTimestep = 1E7
-time.startTimestep = 1E4
-
-adaptivity.useAdaptivity = true
-adaptivity.markingStrategy = threshold
-adaptivity.refinementFraction = 1E-7, 1E-8 | expand prec
-adaptivity.coarseningFraction = 1E-8, 1E-9 | expand prec
-adaptivity.maxLevel = 99
-adaptivity.threshold = 1E-16
-
-grid.FEorder = 1
-grid.gridType = rectangular
-grid.initialLevel = 2
-grid.cells = 1 10
-
-[_ode]
-flux = -3e-6
-head_abstol = 3E-4, 4E-5 | expand prec
-flux_abstol = 4E-8, 1E-8 | expand prec
-
-[NewtonParameters]
-AbsoluteLimit = 1E-10
-Reduction = 1E-10
\ No newline at end of file

testing/ode_heterogeneous_2_adapt.mini.in

-include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
-
-__name = ode_heterogeneous_2_adapt
-_test_command = run
-_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
-_evaluation = ode
-
-output.fileName = ode_heterogeneous_2_adapt | unique
-output.outputPath = ode_heterogeneous_2_adapt | unique
-output.verbose = 0
-
-boundary.file = "{_asset_path}/bcs/infiltration_2d.dat"
-
-parameters.arrayFile = "{_asset_path}/parfields/layered.h5"
-
-time.end = 1E7
-time.maxTimestep = 1E7
-time.startTimestep = 1E4
-
-adaptivity.useAdaptivity = true
-adaptivity.markingStrategy = threshold
-adaptivity.refinementFraction = 1E-7, 1E-8 | expand prec
-adaptivity.coarseningFraction = 1E-8, 1E-9 | expand prec
-adaptivity.maxLevel = 99
-adaptivity.threshold = 1E-16
-
-grid.FEorder = 2
-grid.gridType = rectangular
-grid.initialLevel = 2
-grid.cells = 1 10
-
-[_ode]
-flux = -3e-6
-head_abstol = 3E-5, 5E-6 | expand prec
-flux_abstol = 3E-8, 3E-8 | expand prec
-
-[NewtonParameters]
-AbsoluteLimit = 1E-10
-Reduction = 1E-10
\ No newline at end of file

testing/ode_heterogeneous_3_adapt.mini.in

-include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
-
-__name = ode_heterogeneous_3_adapt
-_test_command = run
-_asset_path = "${CMAKE_CURRENT_LIST_DIR}"