Sucessfully solve heat equation!

python/dune/perftool/pdelab/driver.py

@@ -926,32 +926,30 @@ def name_stationarylinearproblemsolver():
 
 
 @preamble
-def typedef_stationarynonlinearproblemsolver(name):
+def typedef_stationarynonlinearproblemsolver(name, go_type):
     include_file("dune/pdelab/newton/newton.hh", filetag="driver")
-    go_type = type_gridoperator(_driver_data['formdata'])
     ls_type = type_linearsolver()
     x_type = type_vector(_driver_data['formdata'])
     return "typedef Dune::PDELab::Newton<{}, {}, {}> {};".format(go_type, ls_type, x_type, name)
 
 
 @symbol
-def type_stationarynonlinearproblemssolver():
-    typedef_stationarynonlinearproblemsolver("SNP")
+def type_stationarynonlinearproblemssolver(go_type):
+    typedef_stationarynonlinearproblemsolver("SNP", go_type)
     return "SNP"
 
 
 @preamble
-def define_stationarynonlinearproblemsolver(name):
-    snptype = type_stationarynonlinearproblemssolver()
-    go = name_gridoperator(_driver_data['formdata'])
+def define_stationarynonlinearproblemsolver(name, go_type, go):
+    snptype = type_stationarynonlinearproblemssolver(go_type)
     x = name_vector(_driver_data['formdata'])
     ls = name_linearsolver()
     return "{} {}({}, {}, {});".format(snptype, name, go, x, ls)
 
 
 @symbol
-def name_stationarynonlinearproblemsolver():
-    define_stationarynonlinearproblemsolver("snp")
+def name_stationarynonlinearproblemsolver(go_type, go):
+    define_stationarynonlinearproblemsolver("snp", go_type, go)
     return "snp"
 
 
@@ -1015,7 +1013,7 @@ def name_instationarygridoperator():
 def typedef_onestepmethod(name):
     r_type = type_range()
     igo_type = type_instationarygridoperator()
-    snp_type = type_stationarynonlinearproblemssolver()
+    snp_type = type_stationarynonlinearproblemssolver(igo_type)
     vector_type = type_vector(_driver_data['formdata'])
     return "typedef Dune::PDELab::OneStepMethod<{}, {}, {}, {}, {}> {};".format(r_type, igo_type, snp_type, vector_type, vector_type, name)
 
@@ -1031,8 +1029,9 @@ def define_onestepmethod(name):
     ilptype = type_onestepmethod()
     explicit = False
     tsm = name_timesteppingmethod(explicit)
+    igo_type = type_instationarygridoperator()
     igo = name_instationarygridoperator()
-    snp = name_stationarynonlinearproblemsolver()
+    snp = name_stationarynonlinearproblemsolver(igo_type, igo)
     return "{} {}({},{},{});".format(ilptype, name, tsm, igo, snp)
 
 
@@ -1073,36 +1072,6 @@ def name_explicitonestepmethod():
     return "eosm"
 
 
-def time_loop(osm):
-    ini = name_initree()
-    params = name_parameters()
-    formdata = _driver_data['formdata']
-    vector_type = type_vector(formdata)
-    vector_name = name_vector(formdata)
-    expr = _driver_data['form'].coefficients()[0].ufl_element()
-    boundary_name = name_boundary_function(expr)
-    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.01);".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(osm, vector_name, boundary_name, vector_name),
-            "",
-            "  {} = {}new;".format(vector_name, vector_name),
-            "  time += dt;",
-            "}",
-            ""]
-
-
 def explicit_time_loop(eosm):
     ini = name_initree()
     formdata = _driver_data['formdata']
@@ -1116,15 +1085,14 @@ def explicit_time_loop(eosm):
 
     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),
             "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),
+            "  {}.apply(time, dt, {}, {}new);".format(eosm, vector_name, vector_name),
             "",
             "  {} = {}new;".format(vector_name, vector_name),
             "  time += dt;",
@@ -1137,29 +1105,20 @@ def dune_solve():
     # Test if form is linear in ansatzfunction
     linear = is_linear(_driver_data['form'])
 
-    # Test if problem is stationary
-    stationary = is_stationary()
-
     # Test wether we want to do matrix free operator evaluation
     matrix_free = get_option('matrix_free')
 
-    if linear and stationary and matrix_free:
+    if linear and matrix_free:
         go = name_gridoperator(_driver_data['formdata'])
         x = name_vector()
         include_file("dune/perftool/matrixfree.hh", filetag="driver")
         return "solveMatrixFree({},{});".format(go, x)
-    elif linear and stationary and not matrix_free:
+    elif linear and not matrix_free:
         slp = name_stationarylinearproblemsolver()
         return "{}.apply();".format(slp)
-    elif linear and not stationary and not matrix_free:
-        # TODO -> form compiler argument for switching explicit/implicit
-        # osm = name_onestepmethod()
-        # return time_loop(osm)
-        eosm = name_explicitonestepmethod()
-        return explicit_time_loop(eosm)
-    elif not linear and stationary and matrix_free:
+    elif not linear and matrix_free:
         raise NotImplementedError("Nonlinear and matrix free is not yet implemented")
-    elif not linear and stationary and not matrix_free:
+    elif not linear and not matrix_free:
         snp = name_stationarynonlinearproblemsolver()
         return "{}.apply();".format(snp)
     else:
@@ -1328,13 +1287,138 @@ def vtkoutput():
             "{}.write({}, Dune::VTK::ascii);".format(vtkwriter, vtkfile)]
 
 
+
+@preamble
+def define_time(name):
+    return "double {} = 0.0;".format(name)
+
+
+@symbol
+def name_time():
+    define_time("time")
+    return "time"
+
+
+@preamble
+def typedef_vtk_sequence_writer(name):
+    include_file("dune/grid/io/file/vtk/vtksequencewriter.hh", filetag="driver")
+    gv_type = type_leafview()
+    return "typedef Dune::VTKSequenceWriter<{}> {};".format(gv_type, name)
+
+
+@symbol
+def type_vtk_sequence_writer():
+    typedef_vtk_sequence_writer("VTKSW")
+    return "VTKSW"
+
+
+@preamble
+def define_vtk_sequence_writer(name):
+    vtksw_type = type_vtk_sequence_writer()
+    vtkw_type = type_vtkwriter()
+    vtkw = name_vtkwriter()
+    vtkfile = name_vtkfile()
+    return "{} {}(std::make_shared<{}>({}), {});".format(vtksw_type, name, vtkw_type, vtkw, vtkfile)
+
+
+@symbol
+def name_vtk_sequence_writer():
+    define_vtk_sequence_writer("vtkSequenceWriter")
+    return "vtkSequenceWriter"
+
+
+@preamble
+def visualize_initial_condition():
+    include_file("dune/pdelab/gridfunctionspace/vtk.hh", filetag="driver")
+    vtkwriter = name_vtk_sequence_writer()
+    element = _driver_data['form'].coefficients()[0].ufl_element()
+    define_gfs_name(element)
+    gfs = name_gfs(element)
+    vector = name_vector(_driver_data['formdata'])
+    predicate = name_predicate()
+    time = name_time()
+    return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vector, predicate),
+            "{}.write({}, Dune::VTK::appendedraw);".format(vtkwriter, time)]
+
+
+@preamble
+def implicit_time_loop():
+    ini = name_initree()
+    formdata = _driver_data['formdata']
+    params = name_parameters(formdata)
+    time = name_time()
+    expr = _driver_data['form'].coefficients()[0].ufl_element()
+    bctype = name_bctype_function(expr)
+    gfs = name_gfs(expr)
+    cc = name_constraintscontainer(expr)
+    vector_type = type_vector(formdata)
+    vector = name_vector(formdata)
+    osm = name_onestepmethod()
+    boundary = name_boundary_function(expr)
+    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),
+            "while (time<T-1e-8){",
+            "  // Assemble constraints for new time step",
+            "  {}.setTime({}+dt);".format(params, time),
+            "  Dune::PDELab::constraints({}, {}, {});".format(bctype, gfs, cc),
+            "",
+            "  // Do time step",
+            "  {} {}new({});".format(vector_type, vector, vector),
+            "  {}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector),
+            "",
+            "  // Accept new time step",
+            "  {} = {}new;".format(vector, vector),
+            "  time += dt;",
+            "",
+            "  // Output to VTK File",
+            "  {}.write({}, Dune::VTK::appendedraw);".format(vtk_sequence_writer, time),
+            "}",
+            ""]
+
+
+def time_loop():
+    # Test if form is linear in ansatzfunction
+    linear = is_linear(_driver_data['form'])
+
+    # 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:
+        assert False
+    elif linear and not matrix_free and not explicit:
+        implicit_time_loop()
+    if not linear and matrix_free and explicit:
+        assert False
+    elif not linear and matrix_free and not explicit:
+        assert False
+    elif not linear and not matrix_free and explicit:
+        assert False
+    elif not linear and not matrix_free and not explicit:
+        assert False
+
+
 def generate_driver(formdatas, data):
     # The driver module uses a global dictionary for storing necessary data
     set_driver_data(formdatas, data)
 
     # The vtkoutput is the generating method that triggers all others.
     # Alternatively, one could use the `solve` method.
-    vtkoutput()
+    if is_stationary():
+        vtkoutput()
+    else:
+        time_loop()
 
     from dune.perftool.generation import retrieve_cache_items
     from cgen import FunctionDeclaration, FunctionBody, Block, Value

test/heatequation/driver.hh

+/********************************************************/
+// Beware of line number changes, they may corrupt docu!
+//! \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)
+{
+  // dimension and important types
+  const int dim = GV::dimension;
+  typedef double RF;                   // type for computations
+
+  // make user functions and set initial time
+  RF time = 0.0;
+  RF eta = ptree.get("problem.eta",(RF)1.0);
+  Problem<RF> problem(eta);
+  problem.setTime(time);
+  auto glambda = [&](const auto& x){ Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2()); };
+  // auto glambda = [&](const auto& e, const auto& x)
+  //   {return problem.g(e,x);};
+  auto g = Dune::PDELab::
+    makeInstationaryGridFunctionFromCallable(gv,glambda,problem);
+  auto blambda = [&](const auto& i, const auto& x)
+    {return true;};
+  // auto blambda = [&](const auto& i, const auto& x)
+  //   {return problem.b(i,x);};
+  auto b = Dune::PDELab::
+    makeBoundaryConditionFromCallable(gv,blambda);
+
+  // Make grid function space
+  typedef Dune::PDELab::ConformingDirichletConstraints CON;
+  typedef Dune::PDELab::istl::VectorBackend<> VBE;
+  typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE> GFS;
+  GFS gfs(gv,fem);
+  gfs.name("Vh");
+
+  // Assemble constraints
+  typedef typename GFS::template
+    ConstraintsContainer<RF>::Type CC;
+  CC cc;
+  Dune::PDELab::constraints(b,gfs,cc); // assemble constraints
+  std::cout << "constrained dofs=" << cc.size() << " of "
+            << gfs.globalSize() << std::endl;
+
+  // A coefficient vector
+  using Z = Dune::PDELab::Backend::Vector<GFS,RF>;
+  Z z(gfs); // initial value
+
+  // Make a grid function out of it
+  typedef Dune::PDELab::DiscreteGridFunction<GFS,Z> ZDGF;
+  ZDGF zdgf(gfs,z);
+
+  // initialize simulation time,  the coefficient vector
+  Dune::PDELab::interpolate(g,gfs,z);
+
+
+  // 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 Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
+  int degree = ptree.get("fem.degree",(int)1);
+  MBE mbe((int)pow(1+2*degree,dim));
+  typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,
+                                     RF,RF,RF,CC,CC> GO0;
+  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 Dune::PDELab::GridOperator<GFS,GFS,TLOP,MBE,
+                                     RF,RF,RF,CC,CC> GO1;
+  GO1 go1(gfs,cc,gfs,cc,tlop,mbe);
+  typedef Dune::PDELab::OneStepGridOperator<GO0,GO1> IGO;
+  IGO igo(go0,go1);
+
+  // Select a linear solver backend
+  typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<IGO> LS;
+  LS ls(100,0);
+
+  // solve nonlinear problem
+  typedef Dune::PDELab::Newton<IGO,LS,Z> PDESOLVER;
+  PDESOLVER newton(igo,z,ls);
+  newton.setReassembleThreshold(0.0);
+  newton.setVerbosityLevel(2);
+  newton.setReduction(1e-8);
+  newton.setMinLinearReduction(1e-4);
+  newton.setMaxIterations(25);
+  newton.setLineSearchMaxIterations(10);
+
+  // select and prepare time-stepping scheme
+  Dune::PDELab::OneStepThetaParameter<RF> method1(1.0);
+  Dune::PDELab::Alexander2Parameter<RF> method2;
+  Dune::PDELab::Alexander3Parameter<RF> method3;
+  int torder = ptree.get("fem.torder",(int)1);
+  Dune::PDELab::TimeSteppingParameterInterface<RF>*
+    pmethod=&method1;
+  if (torder==1) pmethod = &method1;
+  if (torder==2) pmethod = &method2;
+  if (torder==3) pmethod = &method3;
+  if (torder<1||torder>3)
+    std::cout<<"torder not in [1,3]"<<std::endl;
+  typedef Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVER,Z,Z> OSM;
+  OSM  osm(*pmethod,igo,newton);
+  osm.setVerbosityLevel(2);
+
+  // prepare VTK writer and write first file
+  std::string filename=ptree.get("output.filename","output");
+  struct stat st;
+  if( stat( filename.c_str(), &st ) != 0 )
+    {
+      int stat = 0;
+      stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+      if( stat != 0 && stat != -1)
+        std::cout << "Error: Cannot create directory "
+                  << filename << std::endl;
+    }
+  int subsampling=ptree.get("output.subsampling",(int)0);
+  typedef Dune::SubsamplingVTKWriter<GV> VTKWRITER;
+  VTKWRITER vtkwriter(gv,subsampling);
+  typedef Dune::VTKSequenceWriter<GV> VTKSEQUENCEWRITER;
+  VTKSEQUENCEWRITER vtkSequenceWriter(
+    std::make_shared<VTKWRITER>(vtkwriter),filename,filename,"");
+
+
+  // typedef Dune::PDELab::VTKGridFunctionAdapter<ZDGF> VTKF;
+  // vtkSequenceWriter.addVertexData(std::shared_ptr<VTKF>(
+  //                            new VTKF(zdgf,"solution")));
+  CuttingPredicate predicate;
+  Dune::PDELab::addSolutionToVTKWriter(vtkSequenceWriter, gfs, z, Dune::PDELab::vtk::defaultNameScheme(), predicate);
+  vtkSequenceWriter.write(time,Dune::VTK::appendedraw);
+
+
+  // time loop
+  RF T = ptree.get("problem.T",(RF)1.0);
+  RF dt = ptree.get("fem.dt",(RF)0.1);
+  while (time<T-1e-8)
+    {
+      // assemble constraints for new time step
+      problem.setTime(time+dt);
+      Dune::PDELab::constraints(b,gfs,cc);
+
+      // do time step
+      Z znew(z);
+      osm.apply(time,dt,z,g,znew);
+
+      // accept time step
+      z = znew;
+      time+=dt;
+
+      // output to VTK file
+      vtkSequenceWriter.write(time,Dune::VTK::appendedraw);
+    }
+}