compile Simulation instantiations into library, link to executable

* added 'impl' folder containing explicit template instantiations
* moved Simulation member implementations into .cc file

CMakeLists.txt

@@ -36,6 +36,7 @@ dune_require_cxx_standard(MODULE "dorie" VERSION 14)
 add_subdirectory("bin")
 add_subdirectory("m4")
 add_subdirectory("dune")
+add_subdirectory("lib")
 add_subdirectory("python")
 add_subdirectory("doc")
 if(dune-testtools_FOUND)

dune/dorie/CMakeLists.txt

 add_subdirectory(interface)
 add_subdirectory(solver)
+add_subdirectory(impl)
 
 add_executable("dorie" dorie.cc)
-target_link_dune_default_libraries("dorie")
+#target_link_dune_default_libraries("dorie")
+dune_target_link_libraries(dorie doriesim)
 
 add_executable("sim" EXCLUDE_FROM_ALL sim.cc)
 target_link_dune_default_libraries("sim")
\ No newline at end of file

dune/dorie/dorie.cc

@@ -5,7 +5,6 @@
 #endif
 #include <dune/dorie/dorie.hh> // System and DUNE Headers
 
-
 //===============================================================
 // Main program with grid setup
 //===============================================================
@@ -22,9 +21,13 @@
 * \see Dune::Dorie::FlowSource Class handling Source Term query functions
 */
 
-/// Resolve the second (default) template parameter of YaspGrid
-template<int dim>
-using YaspGrid = Dune::YaspGrid<dim>;
+namespace Dune{
+namespace Dorie{
+	/// Resolve the second (default) template parameter of YaspGrid
+	template<int dim>
+	using YaspGrid = Dune::YaspGrid<dim>;
+}
+}
 
 template<typename Traits>
 using Sim = Dune::Dorie::Simulation<Traits>;
@@ -38,7 +41,7 @@ using SimplexAdaptive = Dune::Dorie::BaseTraits<Dune::UGGrid,
 	Dune::GeometryType::BasicType::simplex,dim,order,true,true>;
 
 template<int dim, int order>
-using Cube = Dune::Dorie::BaseTraits<YaspGrid,
+using Cube = Dune::Dorie::BaseTraits<Dune::Dorie::YaspGrid,
 	Dune::GeometryType::BasicType::cube,dim,order,true,false>;
 
 template<int dim, int order>

dune/dorie/impl/CMakeLists.txt

+dune_add_library(simulations OBJECT 
+	sim_ug_2_1.cc
+	sim_ug_2_2.cc
+	sim_ug_2_3.cc
+	sim_ug_3_1.cc
+	sim_ug_3_2.cc
+	sim_ug_3_3.cc
+	sim_yasp_2_1.cc
+	sim_yasp_2_2.cc
+	sim_yasp_2_3.cc
+	sim_yasp_3_1.cc
+	sim_yasp_3_2.cc
+	sim_yasp_3_3.cc)
\ No newline at end of file

dune/dorie/impl/impl.hh

+#ifndef DUNE_DORIE_IMPL_HH
+#define DUNE_DORIE_IMPL_HH
+
+namespace Dune{
+namespace Dorie{
+
+template<int dim>
+using YaspGrid = Dune::YaspGrid<dim>;
+
+using Geo = Dune::GeometryType::BasicType;
+
+}
+}
+
+#endif // DUNE_DORIE_IMPL_HH

dune/dorie/impl/sim_ug_2_1.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<UGGrid,Geo::cube,2,1,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,2,1,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,2,1,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_ug_2_2.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<UGGrid,Geo::cube,2,2,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,2,2,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,2,2,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_ug_2_3.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<UGGrid,Geo::cube,2,3,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,2,3,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,2,3,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_ug_3_1.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<UGGrid,Geo::cube,3,1,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,3,1,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,3,1,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_ug_3_2.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<UGGrid,Geo::cube,3,2,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,3,2,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,3,2,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_ug_3_3.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<UGGrid,Geo::cube,3,3,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,3,3,true,true>>;
+template class Simulation<BaseTraits<UGGrid,Geo::simplex,3,3,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_yasp_2_1.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<YaspGrid,Geo::cube,2,1,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_yasp_2_2.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<YaspGrid,Geo::cube,2,2,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_yasp_2_3.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<YaspGrid,Geo::cube,2,3,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_yasp_3_1.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<YaspGrid,Geo::cube,3,1,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_yasp_3_2.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<YaspGrid,Geo::cube,3,2,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/impl/sim_yasp_3_3.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dorie/dorie.hh>
+#include <dune/dorie/impl/impl.hh>
+#include <dune/dorie/interface/simulation.cc>
+
+namespace Dune{
+namespace Dorie{
+
+template class Simulation<BaseTraits<YaspGrid,Geo::cube,3,3,true,false>>;
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/interface/simulation.cc

-#include "simulation_traits.hh"
+#include "simulation.hh"
 
-extern template Dune::Dorie::SimulationBase<Dune::Dorie::BaseTraits<Dune::UGGrid,Dune::GeometryType::BasicType::cube,2,1>>;
-extern template Dune::Dorie::SimulationBase<Dune::Dorie::BaseTraits<Dune::UGGrid,Dune::GeometryType::BasicType::cube,2,2>>;
-extern template Dune::Dorie::SimulationBase<Dune::Dorie::BaseTraits<Dune::UGGrid,Dune::GeometryType::BasicType::cube,2,3>>;
\ No newline at end of file
+namespace Dune{
+namespace Dorie{
+
+template<typename Traits>
+Simulation<Traits>::Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile) :
+	helper(_helper), grid(_grid), inifile(_inifile), gv(grid->leafGridView()),
+	mbe_slop(estimate_mbe_entries<RF>(Traits::dim,Traits::GridGeometryType)),
+	mbe_tlop(1.0),
+	verbose(inifile.get<int>("output.verbose"))
+{
+	Dune::Timer timer;
+
+	// --- Grid Function Space ---
+	gfs = std::make_unique<GFS>(GFSHelper::create(gv));
+	gfs->name("matric_head");
+	cc = std::make_unique<CC>();
+	Dune::PDELab::constraints(*gfs,*cc,true);
+
+	// --- Parameter class ---
+	InterpolatorFactory intfac(inifile,helper,verbose);
+	interp = intfac.create();
+	ParameterizationFactory pfac(inifile,helper,*interp,verbose);
+	param = pfac.create();
+
+	// --- Operator Helper Classes ---
+	fboundary = std::make_unique<FlowBoundary>(inifile);
+	fsource = std::make_unique<FlowSource>(inifile,*param);
+	finitial = std::make_unique<FlowInitial>(inifile,*param,gv);
+
+	// --- Local Operators ---
+	slop = std::make_unique<SLOP>(inifile,*param,gv,*fboundary,*fsource);
+	tlop = std::make_unique<TLOP>(inifile,*param,gv);
+	controller = std::make_unique<CalculationController>(inifile,*fboundary,helper);
+
+	// --- Solution Vectors and Initial Condition ---
+	uold = std::make_unique<U>(*gfs,.0);
+	unew = std::make_unique<U>(*gfs,.0);
+	Dune::PDELab::interpolate(*finitial,*gfs,*uold);
+	*unew = *uold;
+
+	// --- Operator Setup --- //
+	operator_setup();
+
+	// --- Utility Class Setup --- //
+	output = OutputWriterFactory::create(inifile,gv,*gfs,*param,*unew);
+
+	AdaptivityHandlerFactory adaptivity_fac(inifile,*grid);
+	adaptivity = adaptivity_fac.create();
+
+	if(verbose>1 && helper.rank()==0)
+			std::cout << ":::  setup time " << std::setw(12) << std::setprecision(4) << std::scientific << timer.elapsed() << std::endl;
+}
+
+template<typename Traits>
+void Simulation<Traits>::operator_setup ()
+{
+	Dune::PDELab::constraints(*this->gfs,*this->cc,false);
+	// --- Grid Operators ---
+	go0 = std::make_unique<GO0>(*gfs,*cc,*gfs,*cc,*slop,mbe_slop);
+	go1 = std::make_unique<GO1>(*gfs,*cc,*gfs,*cc,*tlop,mbe_tlop);
+	igo = std::make_unique<IGO>(*go0,*go1);
+
+	// --- Solvers ---
+	ls = std::make_unique<LS>(5000,0);
+	pdesolver = std::make_unique<PDESOLVER>(*igo,*ls);
+	pdesolver->setParameters(inifile.sub("NewtonParameters"));
+	pdesolver->setVerbosityLevel(verbose);
+
+	// --- Time Step Operators ---
+	Dune::Dorie::OneStepMethodResult osm_result;
+	if(osm){
+		osm_result = osm->result(); // cache old result
+	}
+	osm = std::make_unique<OSM>(alex2,*igo,*pdesolver,osm_result);
+	osm->setVerbosityLevel(verbose+1);
+	osm->setStepNumber(osm_result.successful.timesteps+1);
+
+	gfs->update();
+
+	// some output makes it look intelligent
+	if(verbose>0){
+		ios_base_all_saver restore(std::cout);
+		if (helper.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 " << helper.rank() << ": number of DOF: " << DOFnumber << std::endl;
+		DOFnumber = gv.comm().sum(DOFnumber);
+		if (helper.rank()==0)
+			std::cout << "  Total number of DOF: " << DOFnumber << std::endl;
+	}
+}
+
+template<typename Traits>
+bool Simulation<Traits>::compute_time_step ()
+{
+	const RF t = controller->getTime();
+	const RF dt = controller->getDT();
+	pdesolver->setMaxIterations(controller->getIterations());
+	bool exception = false;
+	try
+	{
+		osm->apply(t,dt,*uold,*unew);
+		*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 reacts to outcome of solution
+	return controller->validate(exception);
+}
+
+template<typename Traits>
+void Simulation<Traits>::run ()
+{
+	output->write_vtk_output(controller->getTime());
+	while(controller->doStep()) {
+		if(!compute_time_step()){
+			continue;
+		}
+		if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *uold, *unew)){ // reset operators if grid changes
+			operator_setup();
+		}
+		output->write_vtk_output(controller->getTime());
+	}
+}
+
+} // namespace Dorie
+} // namespace Dune
\ No newline at end of file

dune/dorie/interface/simulation.hh

-#ifndef DUNE_DORIE_SIMULATION_BASE_HH
-#define DUNE_DORIE_SIMULATION_BASE_HH
+#ifndef DUNE_DORIE_SIMULATION_HH
+#define DUNE_DORIE_SIMULATION_HH
 
 namespace Dune{
 namespace Dorie{
@@ -112,19 +112,7 @@ protected:
 public:
 
 	/// Execute the simulation until tEnd is reached.
-	virtual void run ()
-	{
-		output->write_vtk_output(controller->getTime());
-		while(controller->doStep()) {
-			if(!compute_time_step()){
-				continue;
-			}
-			if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *uold, *unew)){ // reset operators if grid changes
-				operator_setup();
-			}
-			output->write_vtk_output(controller->getTime());
-		}
-	}
+	void run ();
 
 	/// Construct basic simulation. Prepare setup for basic simulation.
 	/**
@@ -132,136 +120,23 @@ public:
 	 *  \param _grid Shared pointer to the grid
 	 *  \param _inifile Dune parameter file parser
 	 */
-	 Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile) :
-		helper(_helper), grid(_grid), inifile(_inifile), gv(grid->leafGridView()),
-		mbe_slop(estimate_mbe_entries<RF>(Traits::dim,Traits::GridGeometryType)),
-		mbe_tlop(1.0),
-		verbose(inifile.get<int>("output.verbose"))
-	{
-		Dune::Timer timer;
-
-		// --- Grid Function Space ---
-		gfs = std::make_unique<GFS>(GFSHelper::create(gv));
-		gfs->name("matric_head");
-		cc = std::make_unique<CC>();
-		Dune::PDELab::constraints(*gfs,*cc,true);
-
-		// --- Parameter class ---
-		InterpolatorFactory intfac(inifile,helper,verbose);
-		interp = intfac.create();
-		ParameterizationFactory pfac(inifile,helper,*interp,verbose);
-		param = pfac.create();
-
-		// --- Operator Helper Classes ---
-		fboundary = std::make_unique<FlowBoundary>(inifile);
-		fsource = std::make_unique<FlowSource>(inifile,*param);
-		finitial = std::make_unique<FlowInitial>(inifile,*param,gv);
-
-		// --- Local Operators ---
-		slop = std::make_unique<SLOP>(inifile,*param,gv,*fboundary,*fsource);
-		tlop = std::make_unique<TLOP>(inifile,*param,gv);
-		controller = std::make_unique<CalculationController>(inifile,*fboundary,helper);
-
-		// --- Solution Vectors and Initial Condition ---
-		uold = std::make_unique<U>(*gfs,.0);
-		unew = std::make_unique<U>(*gfs,.0);
-		Dune::PDELab::interpolate(*finitial,*gfs,*uold);
-		*unew = *uold;
-
-		// --- Operator Setup --- //
-		operator_setup();
-
-		// --- Utility Class Setup --- //
-		output = OutputWriterFactory::create(inifile,gv,*gfs,*param,*unew);
-
-		AdaptivityHandlerFactory adaptivity_fac(inifile,*grid);
-		adaptivity = adaptivity_fac.create();
-
-		if(verbose>1 && helper.rank()==0)
-				std::cout << ":::  setup time " << std::setw(12) << std::setprecision(4) << std::scientific << timer.elapsed() << std::endl;
-	}
+	Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile);
 
 	/// Create the Grid Operators, Solvers, and Time Step Operator
 	/** This function must be called after changes to grid or GFS!
 	 */
-	void operator_setup ()
-	{
-		Dune::PDELab::constraints(*this->gfs,*this->cc,false);
-		// --- Grid Operators ---
-		go0 = std::make_unique<GO0>(*gfs,*cc,*gfs,*cc,*slop,mbe_slop);
-		go1 = std::make_unique<GO1>(*gfs,*cc,*gfs,*cc,*tlop,mbe_tlop);
-		igo = std::make_unique<IGO>(*go0,*go1);
-
-		// --- Solvers ---
-		ls = std::make_unique<LS>(5000,0);
-		pdesolver = std::make_unique<PDESOLVER>(*igo,*ls);
-		pdesolver->setParameters(inifile.sub("NewtonParameters"));
-		pdesolver->setVerbosityLevel(verbose);
-
-		// --- Time Step Operators ---
-		Dune::Dorie::OneStepMethodResult osm_result;
-		if(osm){
-			osm_result = osm->result(); // cache old result
-		}
-		osm = std::make_unique<OSM>(alex2,*igo,*pdesolver,osm_result);
-		osm->setVerbosityLevel(verbose+1);
-		osm->setStepNumber(osm_result.successful.timesteps+1);
-
-		gfs->update();
-
-		// some output makes it look intelligent
-		if(verbose>0){
-			ios_base_all_saver restore(std::cout);
-			if (helper.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 " << helper.rank() << ": number of DOF: " << DOFnumber << std::endl;
-			DOFnumber = gv.comm().sum(DOFnumber);
-			if (helper.rank()==0)
-				std::cout << "  Total number of DOF: " << DOFnumber << std::endl;
-		}
-	}
+	void operator_setup ();
 
 	/// Compute a time step
 	/** Catch errors in the Newton and the Linear Solver.
 	 *  \throw Dune::Exception Fatal error occurs during computation
 	 *  \return True if computation succeeded
 	 */
-	bool compute_time_step ()
-	{
-		const RF t = controller->getTime();
-		const RF dt = controller->getDT();
-		pdesolver->setMaxIterations(controller->getIterations());
-		bool exception = false;
-		try
-		{
-			osm->apply(t,dt,*uold,*unew);
-			*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!" <<