Add explicit onestep methods

python/dune/perftool/options.py

@@ -14,6 +14,7 @@ def get_form_compiler_arguments():
     parser.add_argument('--version', action='version', version='%(prog)s 0.1')
     parser.add_argument("--numerical-jacobian", action="store_true", help="use numerical jacobians (only makes sense, if uflpdelab for some reason fails to generate analytic jacobians)")
     parser.add_argument("--matrix-free", action="store_true", help="use iterative solver with matrix free jacobian application")
+    parser.add_argument("--explicit-time-stepping", action="store_true", help="use explicit time stepping")
     parser.add_argument(
         "--exact-solution-expression",
         type=str,

python/dune/perftool/pdelab/driver.py

@@ -957,7 +957,11 @@ def name_stationarynonlinearproblemsolver(go_type, go):
 @preamble
 def typedef_timesteppingmethod(name):
     r_type = type_range()
-    return "typedef Dune::PDELab::OneStepThetaParameter<{}> {};".format(r_type, name)
+    explicit = get_option('explicit_time_stepping')
+    if explicit:
+        return "typedef Dune::PDELab::ExplicitEulerParameter<{}> {};".format(r_type, name)
+    else:
+        return "typedef Dune::PDELab::OneStepThetaParameter<{}> {};".format(r_type, name)
 
 
 @symbol
@@ -967,18 +971,18 @@ def type_timesteppingmethod():
 
 
 @preamble
-def define_timesteppingmethod(name, explicit):
+def define_timesteppingmethod(name):
     tsm_type = type_timesteppingmethod()
-    # TODO enable theta from ini file
+    explicit = get_option('explicit_time_stepping')
     if explicit:
-        return "{} {}(0.0);".format(tsm_type, name)
+        return "{} {};".format(tsm_type, name)
     else:
         return "{} {}(1.0);".format(tsm_type, name)
 
 
 @symbol
-def name_timesteppingmethod(explicit):
-    define_timesteppingmethod("tsm", explicit)
+def name_timesteppingmethod():
+    define_timesteppingmethod("tsm")
     return "tsm"
 
 
@@ -987,7 +991,11 @@ def typedef_instationarygridoperator(name):
     include_file("dune/pdelab/gridoperator/onestep.hh", filetag="driver")
     go_type = type_gridoperator(_driver_data['formdata'])
     mass_go_type = type_gridoperator(_driver_data['mass_formdata'])
-    return "typedef Dune::PDELab::OneStepGridOperator<{},{}> {};".format(go_type, mass_go_type, name)
+    explicit = get_option('explicit_time_stepping')
+    if explicit:
+        return "typedef Dune::PDELab::OneStepGridOperator<{},{},false> {};".format(go_type, mass_go_type, name)
+    else:
+        return "typedef Dune::PDELab::OneStepGridOperator<{},{}> {};".format(go_type, mass_go_type, name)
 
 
 @symbol
@@ -1028,8 +1036,7 @@ def type_onestepmethod():
 @preamble
 def define_onestepmethod(name):
     ilptype = type_onestepmethod()
-    explicit = False
-    tsm = name_timesteppingmethod(explicit)
+    tsm = name_timesteppingmethod()
     igo_type = type_instationarygridoperator()
     igo = name_instationarygridoperator()
     snp = name_stationarynonlinearproblemsolver(igo_type, igo)
@@ -1060,8 +1067,7 @@ def type_explicitonestepmethod():
 @preamble
 def define_explicitonestepmethod(name):
     eosm_type = type_explicitonestepmethod()
-    explicit = True
-    tsm = name_timesteppingmethod(explicit)
+    tsm = name_timesteppingmethod()
     igo = name_instationarygridoperator()
     ls = name_linearsolver()
     return "{} {}({}, {}, {});".format(eosm_type, name, tsm, igo, ls)
@@ -1073,34 +1079,6 @@ def name_explicitonestepmethod():
     return "eosm"
 
 
-def explicit_time_loop(eosm):
-    ini = name_initree()
-    formdata = _driver_data['formdata']
-    params = name_parameters(formdata)
-    vector_type = type_vector(formdata)
-    vector_name = name_vector(formdata)
-    expr = _driver_data['form'].coefficients()[0].ufl_element()
-    bctype_name = name_bctype_function(expr)
-    gfs_name = name_gfs(expr)
-    cc_name = name_constraintscontainer(expr)
-
-    return ["",
-            "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
-            "double dt = {}.get<double>(\"instat.dt\", 0.1);".format(ini),
-            "double time = 0;",
-            "while (time<T-1e-8){",
-            "  {}.setTime(time+dt);".format(params),
-            "  Dune::PDELab::constraints({}, {}, {});".format(bctype_name, gfs_name, cc_name),
-            "",
-            "  {} {}new({});".format(vector_type, vector_name, vector_name),
-            "  {}.apply(time, dt, {}, {}new);".format(eosm, vector_name, vector_name),
-            "",
-            "  {} = {}new;".format(vector_name, vector_name),
-            "  time += dt;",
-            "}",
-            ""]
-
-
 @preamble
 def dune_solve():
     # Test if form is linear in ansatzfunction
@@ -1357,14 +1335,22 @@ def implicit_time_loop():
     cc = name_constraintscontainer(expr)
     vector_type = type_vector(formdata)
     vector = name_vector(formdata)
-    osm = name_onestepmethod()
-    boundary = name_boundary_function(expr)
+
+    # TODO -> formcompiler argument
+    explicit = get_option('explicit_time_stepping')
+    if explicit:
+        osm = name_explicitonestepmethod()
+        apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
+    else:
+        boundary = name_boundary_function(expr)
+        osm = name_onestepmethod()
+        apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
     visualize_initial_condition()
     vtk_sequence_writer = name_vtk_sequence_writer()
 
     return ["",
             "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
-            "double dt = {}.get<double>(\"instat.dt\", 0.01);".format(ini),
+            "double dt = {}.get<double>(\"instat.dt\", 0.1);".format(ini),
             "while (time<T-1e-8){",
             "  // Assemble constraints for new time step",
             "  {}.setTime({}+dt);".format(params, time),
@@ -1372,7 +1358,7 @@ def implicit_time_loop():
             "",
             "  // Do time step",
             "  {} {}new({});".format(vector_type, vector, vector),
-            "  {}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector),
+            "  {}".format(apply_call),
             "",
             "  // Accept new time step",
             "  {} = {}new;".format(vector, vector),
@@ -1391,25 +1377,14 @@ def time_loop():
     # Test wether we want to do matrix free operator evaluation
     matrix_free = get_option('matrix_free')
 
-    # Test if we use explicit time stepping
-    # TODO
-    explicit = False
-
-    if linear and matrix_free and explicit:
-        assert False
-    elif linear and matrix_free and not explicit:
-        assert False
-    elif linear and not matrix_free and explicit:
+    # Create time loop
+    if linear and matrix_free:
         assert False
-    elif linear and not matrix_free and not explicit:
+    elif linear and not matrix_free:
         implicit_time_loop()
-    if not linear and matrix_free and explicit:
-        assert False
-    elif not linear and matrix_free and not explicit:
+    if not linear and matrix_free:
         assert False
-    elif not linear and not matrix_free and explicit:
-        assert False
-    elif not linear and not matrix_free and not explicit:
+    elif not linear and not matrix_free:
         assert False
 
 

test/heatequation/CMakeLists.txt

 add_generated_executable(UFLFILE heatequation.ufl
-                         TARGET heatequation
+                         TARGET heatequation_implicit
                          FORM_COMPILER_ARGS
                          )
 
 
-dune_add_system_test(TARGET heatequation
-                     INIFILE heatequation.mini
+dune_add_system_test(TARGET heatequation_implicit
+                     INIFILE heatequation_implicit.mini
                      )
 
-add_executable(tutorial03 tutorial03.cc)
-dune_symlink_to_source_files(FILES tutorial03.ini)
+add_generated_executable(UFLFILE heatequation.ufl
+                         TARGET heatequation_explicit
+                         FORM_COMPILER_ARGS --explicit-time-stepping
+                         )
+
+
+dune_add_system_test(TARGET heatequation_explicit
+                     INIFILE heatequation_explicit.mini
+                     )
+
+add_executable(heatequation_implicit_ref heatequation_implicit_ref.cc)
+dune_symlink_to_source_files(FILES heatequation_implicit_ref.ini)
+set_target_properties(heatequation_implicit_ref PROPERTIES EXCLUDE_FROM_ALL 1)
+
+add_executable(heatequation_explicit_ref heatequation_explicit_ref.cc)
+set_target_properties(heatequation_explicit_ref PROPERTIES EXCLUDE_FROM_ALL 1)
 
 # # Add the reference code with the PDELab localoperator that produced
 # # the reference vtk file

test/heatequation/driver.hh

@@ -3,7 +3,6 @@
 //! \brief Driver function to set up and solve the problem
 /********************************************************/
 
-#include </home/rhess/code-gen/dune/build/dune-perftool/test/heatequation/heatequation_operator.hh>
 #include <dune/perftool/vtkpredicate.hh>
 template<typename GV, typename FEM>
 void driver (const GV& gv, const FEM& fem, Dune::ParameterTree& ptree)
@@ -57,11 +56,11 @@ void driver (const GV& gv, const FEM& fem, Dune::ParameterTree& ptree)
 
 
   // Make instationary grid operator
-  // typedef NonlinearHeatFEM<Problem<RF>,FEM> LOP;
-  // LOP lop(problem);
-  using LOP = LocalOperatorPoisson<GFS, GFS>;
-  LocalOperatorPoissonParams params_poisson;
-  LOP lop(ptree, params_poisson);
+  typedef NonlinearHeatFEM<Problem<RF>,FEM> LOP;
+  LOP lop(problem);
+  // using LOP = LocalOperatorPoisson<GFS, GFS>;
+  // LocalOperatorPoissonParams params_poisson;
+  // LOP lop(ptree, params_poisson);
 
 
   typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
@@ -72,11 +71,11 @@ void driver (const GV& gv, const FEM& fem, Dune::ParameterTree& ptree)
   GO0 go0(gfs,cc,gfs,cc,lop,mbe);
 
 
-  // typedef L2<FEM> TLOP;
-  // TLOP tlop;
-  using TLOP = LocalOperatorMass<GFS, GFS>;
-  LocalOperatorMassParams params_mass;
-  TLOP tlop(ptree, params_mass);
+  typedef L2<FEM> TLOP;
+  TLOP tlop;
+  // using TLOP = LocalOperatorMass<GFS, GFS>;
+  // LocalOperatorMassParams params_mass;
+  // TLOP tlop(ptree, params_mass);
 
 
   typedef Dune::PDELab::GridOperator<GFS,GFS,TLOP,MBE,

test/heatequation/heatequation_explicit.mini

+__name = heatequation_explicit
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 8 8
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = heatequation_ref
+extension = vtu
+
+[instat]
+T = 1.0
+dt = 0.01
\ No newline at end of file

test/heatequation/heatequation_explicit_ref.cc

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/pdelab/boilerplate/pdelab.hh>
+#include <dune/pdelab/localoperator/convectiondiffusionfem.hh>
+#include <dune/pdelab/localoperator/l2.hh>
+
+//***********************************************************************
+//***********************************************************************
+// diffusion problem with time dependent coefficients
+//***********************************************************************
+//***********************************************************************
+
+const double kx = 2.0, ky = 2.0;
+
+template<typename GV, typename RF>
+class GenericProblem
+{
+  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
+
+public:
+  typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
+
+  GenericProblem () : time(0.0) {}
+
+  //! tensor diffusion coefficient
+  typename Traits::PermTensorType
+  A (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    typename Traits::PermTensorType I;
+    for (std::size_t i=0; i<Traits::dimDomain; i++)
+      for (std::size_t j=0; j<Traits::dimDomain; j++)
+        I[i][j] = (i==j) ? 1.0 : 0.0;
+    return I;
+  }
+
+  //! velocity field
+  typename Traits::RangeType
+  b (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    typename Traits::RangeType v(0.0);
+    return v;
+  }
+
+  //! sink term
+  typename Traits::RangeFieldType
+  c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! source term
+  typename Traits::RangeFieldType
+  f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! boundary condition type function
+  BCType
+  bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
+  }
+
+  //! Dirichlet boundary condition value
+  typename Traits::RangeFieldType
+  g (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
+  {
+    typename Traits::DomainType x = e.geometry().global(xlocal);
+    Dune::FieldVector<double,2> c(0.5);
+    c-= x;
+    return std::exp(-1.*c.two_norm2());
+    // return std::exp(-(kx*kx+ky*ky)*M_PI*M_PI*time) * sin(kx*M_PI*x[0]) * sin(ky*M_PI*x[1]);
+  }
+
+  //! Neumann boundary condition
+  typename Traits::RangeFieldType
+  j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! Neumann boundary condition
+  typename Traits::RangeFieldType
+  o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! set time for subsequent evaluation
+  void setTime (RF t)
+  {
+    time = t;
+    //std::cout << "setting time to " << time << std::endl;
+  }
+
+private:
+  RF time;
+};
+
+//***********************************************************************
+//***********************************************************************
+// a function that does the simulation on a given grid
+//***********************************************************************
+//***********************************************************************
+
+template<typename GM, unsigned int degree, Dune::GeometryType::BasicType elemtype,
+         Dune::PDELab::MeshType meshtype, Dune::SolverCategory::Category solvertype>
+void do_simulation (double T, double dt, GM& grid, std::string basename)
+{
+  // define parameters
+  typedef double NumberType;
+
+  // make problem parameters
+  typedef GenericProblem<typename GM::LeafGridView,NumberType> Problem;
+  Problem problem;
+  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditionAdapter<Problem> BCType;
+  BCType bctype(grid.leafGridView(),problem);
+
+  // make a finite element space
+  typedef Dune::PDELab::CGSpace<GM,NumberType,degree,BCType,elemtype,meshtype,solvertype> FS;
+  FS fs(grid,bctype);
+
+  // assemblers for finite element problem
+  typedef Dune::PDELab::ConvectionDiffusionFEM<Problem,typename FS::FEM> LOP;
+  LOP lop(problem,4);
+  typedef Dune::PDELab::GalerkinGlobalAssembler<FS,LOP,solvertype> SASS;
+  SASS sass(fs,lop,6);
+  typedef Dune::PDELab::L2 MLOP;
+  MLOP mlop(2*degree);
+  typedef Dune::PDELab::GalerkinGlobalAssembler<FS,MLOP,solvertype> TASS;
+  TASS tass(fs,mlop,6);
+  typedef Dune::PDELab::OneStepGlobalAssembler<SASS,TASS,false> ASSEMBLER;
+  ASSEMBLER assembler(sass,tass);
+
+  // make a degree of freedom vector and set initial value
+  typedef typename FS::DOF V;
+  V x(fs.getGFS(),0.0);
+  typedef Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem> G;
+  G g(grid.leafGridView(),problem);
+  problem.setTime(0.0);
+  Dune::PDELab::interpolate(g,fs.getGFS(),x);
+
+  // linear solver backend
+  typedef Dune::PDELab::ISTLSolverBackend_CG_AMG_SSOR<FS,ASSEMBLER,solvertype> SBE;
+  SBE sbe(fs,assembler,5000,1);
+
+  // linear problem solver
+  typedef Dune::PDELab::StationaryLinearProblemSolver<typename ASSEMBLER::GO,typename SBE::LS,V> PDESOLVER;
+  PDESOLVER pdesolver(*assembler,*sbe,1e-6);
+
+  // // time-stepper
+  // Dune::PDELab::OneStepThetaParameter<NumberType> method(1.0);
+  // typedef Dune::PDELab::OneStepMethod<NumberType,typename ASSEMBLER::GO,PDESOLVER,V> OSM;
+  // OSM osm(method,*assembler,pdesolver);
+  // osm.setVerbosityLevel(2);
+
+  Dune::PDELab::ExplicitEulerParameter<NumberType> method;
+  typedef Dune::PDELab::SimpleTimeController<NumberType> TC;
+  TC tc;
+  typedef Dune::PDELab::ExplicitOneStepMethod<NumberType,typename ASSEMBLER::GO,typename SBE::LS,V,V,TC> OSM;
+  OSM osm(method,*assembler,*sbe,tc);
+  osm.setVerbosityLevel(2);
+
+  // graphics for initial guess
+  Dune::PDELab::FilenameHelper fn(basename);
+  { // start a new block to automatically delete the VTKWriter object
+    Dune::SubsamplingVTKWriter<typename GM::LeafGridView> vtkwriter(grid.leafGridView(),degree-1);
+    typename FS::DGF xdgf(fs.getGFS(),x);
+    vtkwriter.addVertexData(std::make_shared<typename FS::VTKF>(xdgf,"x_h"));
+    vtkwriter.write(fn.getName(),Dune::VTK::appendedraw);
+    fn.increment();
+  }
+
+  // time loop
+  NumberType time = 0.0;
+  while (time<T-1e-10)
+    {
+      // assemble constraints for new time step (assumed to be constant for all substeps)
+      problem.setTime(time+dt);
+      fs.assembleConstraints(bctype);
+
+      // do time step
+      V xnew(fs.getGFS(),0.0);
+      osm.apply(time,dt,x,xnew);
+
+      // output to VTK file
+      {
+        Dune::SubsamplingVTKWriter<typename GM::LeafGridView> vtkwriter(grid.leafGridView(),degree-1);
+        typename FS::DGF xdgf(fs.getGFS(),xnew);
+        vtkwriter.addVertexData(std::make_shared<typename FS::VTKF>(xdgf,"x_h"));
+        vtkwriter.write(fn.getName(),Dune::VTK::appendedraw);
+        fn.increment();
+      }
+
+      // accept time step
+      x = xnew;
+      time += dt;
+    }
+}
+
+//***********************************************************************
+//***********************************************************************
+// the main function
+//***********************************************************************
+//***********************************************************************
+
+int main(int argc, char **argv)
+{
+  // initialize MPI, finalize is done automatically on exit
+  Dune::MPIHelper::instance(argc,argv);
+
+  // read command line arguments
+  if (argc!=4)
+    {
+      std::cout << "usage: " << argv[0] << " <T> <dt> <cells>" << std::endl;
+      return 0;
+    }
+  double T; sscanf(argv[1],"%lg",&T);
+  double dt; sscanf(argv[2],"%lg",&dt);
+  int cells; sscanf(argv[3],"%d",&cells);
+
+  // start try/catch block to get error messages from dune
+  try {
+
+    const int dim=2;
+    const int degree=1;
+    const Dune::SolverCategory::Category solvertype = Dune::SolverCategory::sequential;
+    const Dune::GeometryType::BasicType elemtype = Dune::GeometryType::cube;
+    const Dune::PDELab::MeshType meshtype = Dune::PDELab::MeshType::conforming;
+
+    typedef Dune::YaspGrid<dim> GM;
+    typedef Dune::PDELab::StructuredGrid<GM> Grid;
+    Grid grid(elemtype,cells);
+    grid->loadBalance();
+
+    std::stringstream basename;
+    basename << "heatequation_explicit_ref";
+    do_simulation<GM,degree,elemtype,meshtype,solvertype>(T,dt,*grid,basename.str());
+  }
+  catch (std::exception & e) {
+    std::cout << "STL ERROR: " << e.what() << std::endl;
+    return 1;
+  }
+  catch (Dune::Exception & e) {
+    std::cout << "DUNE ERROR: " << e.what() << std::endl;
+    return 1;
+  }
+  catch (...) {
+    std::cout << "Unknown ERROR" << std::endl;
+    return 1;
+  }
+
+  // done
+  return 0;
+}

test/heatequation/heatequation.mini → test/heatequation/heatequation_implicit.mini

-__name = heatequation
+__name = heatequation_implicit
 
 lowerleft = 0.0 0.0
 upperright = 1.0 1.0

test/heatequation/tutorial03.cc → test/heatequation/heatequation_implicit_ref.cc

@@ -82,7 +82,7 @@ int main(int argc, char** argv)
     // open ini file
     Dune::ParameterTree ptree;
     Dune::ParameterTreeParser ptreeparser;
-    ptreeparser.readINITree("tutorial03.ini",ptree);
+    ptreeparser.readINITree("heatequation_implicit_ref.ini",ptree);
     ptreeparser.readOptions(argc,argv,ptree);
 
     // read ini file

test/heatequation/tutorial03.ini → test/heatequation/heatequation_implicit_ref.ini

 [grid]
 dim=2
-refinement=1
+refinement=3
 
 [grid.structured]
 LX=1.0