From a05252a9c70a1d9b09e972c9b723b504505ad463 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de> Date: Thu, 25 Aug 2016 11:04:34 +0200 Subject: [PATCH] Add explicit onestep methods --- python/dune/perftool/options.py | 1 + python/dune/perftool/pdelab/driver.py | 93 +++---- test/heatequation/CMakeLists.txt | 24 +- test/heatequation/driver.hh | 21 +- test/heatequation/heatequation_explicit.mini | 15 + .../heatequation/heatequation_explicit_ref.cc | 259 ++++++++++++++++++ ...uation.mini => heatequation_implicit.mini} | 2 +- ...rial03.cc => heatequation_implicit_ref.cc} | 2 +- ...al03.ini => heatequation_implicit_ref.ini} | 2 +- 9 files changed, 341 insertions(+), 78 deletions(-) create mode 100644 test/heatequation/heatequation_explicit.mini create mode 100644 test/heatequation/heatequation_explicit_ref.cc rename test/heatequation/{heatequation.mini => heatequation_implicit.mini} (84%) rename test/heatequation/{tutorial03.cc => heatequation_implicit_ref.cc} (98%) rename test/heatequation/{tutorial03.ini => heatequation_implicit_ref.ini} (93%) diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py index bbea79cf..d6761605 100644 --- a/python/dune/perftool/options.py +++ b/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, diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py index 1f8988ae..a3d7e88a 100644 --- a/python/dune/perftool/pdelab/driver.py +++ b/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 diff --git a/test/heatequation/CMakeLists.txt b/test/heatequation/CMakeLists.txt index e0e338b3..3ac71527 100644 --- a/test/heatequation/CMakeLists.txt +++ b/test/heatequation/CMakeLists.txt @@ -1,15 +1,29 @@ 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 diff --git a/test/heatequation/driver.hh b/test/heatequation/driver.hh index e42f73d8..1aaf9f17 100644 --- a/test/heatequation/driver.hh +++ b/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, diff --git a/test/heatequation/heatequation_explicit.mini b/test/heatequation/heatequation_explicit.mini new file mode 100644 index 00000000..2360fb6b --- /dev/null +++ b/test/heatequation/heatequation_explicit.mini @@ -0,0 +1,15 @@ +__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 diff --git a/test/heatequation/heatequation_explicit_ref.cc b/test/heatequation/heatequation_explicit_ref.cc new file mode 100644 index 00000000..28470723 --- /dev/null +++ b/test/heatequation/heatequation_explicit_ref.cc @@ -0,0 +1,259 @@ +// -*- 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; +} diff --git a/test/heatequation/heatequation.mini b/test/heatequation/heatequation_implicit.mini similarity index 84% rename from test/heatequation/heatequation.mini rename to test/heatequation/heatequation_implicit.mini index 00a89305..c49499f7 100644 --- a/test/heatequation/heatequation.mini +++ b/test/heatequation/heatequation_implicit.mini @@ -1,4 +1,4 @@ -__name = heatequation +__name = heatequation_implicit lowerleft = 0.0 0.0 upperright = 1.0 1.0 diff --git a/test/heatequation/tutorial03.cc b/test/heatequation/heatequation_implicit_ref.cc similarity index 98% rename from test/heatequation/tutorial03.cc rename to test/heatequation/heatequation_implicit_ref.cc index c1164ce2..8a486b85 100644 --- a/test/heatequation/tutorial03.cc +++ b/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 diff --git a/test/heatequation/tutorial03.ini b/test/heatequation/heatequation_implicit_ref.ini similarity index 93% rename from test/heatequation/tutorial03.ini rename to test/heatequation/heatequation_implicit_ref.ini index 7844d548..7783f5b0 100644 --- a/test/heatequation/tutorial03.ini +++ b/test/heatequation/heatequation_implicit_ref.ini @@ -1,6 +1,6 @@ [grid] dim=2 -refinement=1 +refinement=3 [grid.structured] LX=1.0 -- GitLab