Commit 6f9304f0 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'update/dune-v2.5.1' into feature/dune-randomfield

parents 9c168d22 8206188f
image: dorie/dune-env:2.5
image: dorie/dune-env:2.5.1
variables:
DUNE_CONTROL_PATH: /opt/dune:$CI_PROJECT_DIR
......@@ -14,7 +14,7 @@ build:main:
tags:
- docker
script:
- MAKE_FLAGS="-j 2" ./dune-common/bin/dunecontrol --only=dorie all
- CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" MAKE_FLAGS="-j 2" ./dune-common/bin/dunecontrol --only=dorie all
- export PATH=/opt/dune/dorie/build-cmake/bin:$PATH
- ARGS="--output-on-failure -j 2" ./dune-common/bin/dunecontrol --only=dorie make test
- ./dune-common/bin/dunecontrol --only=dorie make doc
......@@ -25,25 +25,13 @@ build:main:
- $CI_PROJECT_DIR/build-cmake/doc/html/
expire_in: 1 week
build:python3:
tags:
- docker
script:
- ./dune-common/bin/dunecontrol --only=dune-python exec "rm -rf build-cmake"
- ./dune-common/bin/dunecontrol --only=dune-testtools exec "rm -rf build-cmake"
- ./dune-common/bin/dunecontrol --only=dorie exec "rm -rf build-cmake"
- CMAKE_FLAGS="-DDUNE_FORCE_PYTHON3=True" MAKE_FLAGS="-j 2" ./dune-common/bin/dunecontrol --module=dorie all
- export PATH=/opt/dune/dorie/build-cmake/bin:$PATH
- ARGS="--output-on-failure -j 2" ./dune-common/bin/dunecontrol --only=dorie make test
stage: build
build:update_dune:
tags:
- docker
script:
- ./dune-common/bin/dunecontrol update || true
- ./dune-common/bin/dunecontrol exec "rm -rf build-cmake"
- MAKE_FLAGS="-j 2" ./dune-common/bin/dunecontrol all
- CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" MAKE_FLAGS="-j 2" ./dune-common/bin/dunecontrol all
- export PATH=/opt/dune/dorie/build-cmake/bin:$PATH
- ARGS="--output-on-failure -j 2" ./dune-common/bin/dunecontrol --only=dorie make test
stage: build
......@@ -52,7 +40,7 @@ build:debug:
tags:
- docker
script:
- CMAKE_FLAGS="-DCMAKE_BUILD_TYPE=Debug" ./dune-common/bin/dunecontrol --only=dorie all
- CMAKE_FLAGS="-DCMAKE_BUILD_TYPE=Debug -DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" ./dune-common/bin/dunecontrol --only=dorie all
stage: build
deploy:docs:
......
FROM dorie/dune-env:2.5
FROM dorie/dune-env:2.5.1
MAINTAINER Dion Häfner
WORKDIR /opt/dune
ADD . /opt/dune/dorie/
RUN source /opt/dune/venv/dorie/bin/activate && ./dune-common/bin/dunecontrol --only=dorie all
RUN CMAKE_FLAGS="-DCMAKE_BUILD_TYPE=Release -DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" ./dune-common/bin/dunecontrol --only=dorie all
ENV PATH="/opt/dune/dorie/build-cmake/bin:${PATH}"
WORKDIR /sim
ENTRYPOINT ["/opt/dune/dorie/build-cmake/bin/dorie"]
\ No newline at end of file
This diff is collapsed.
#!/usr/bin/env python
#!/usr/bin/env python3
import os
import sys
......@@ -24,7 +24,7 @@ MPIEXEC_POSTFLAGS = "@MPIEXEC_POSTFLAGS@"
# set some paths
DORIE_EXECUTABLE = os.path.join(DORIEDIR, "dune/dorie/dorie")
PARAMETERDIR = os.path.join(DORIEDIR, "doc/default_files")
DORIE_PYTHON = os.path.join(DORIEDIR, "dune-env")
DORIE_PYTHON = os.path.join(DORIEDIR, "run-in-dune-env")
def MPIRUN(nproc,exe,*args,**kwargs):
mpi_flags = kwargs.get("mpi_flags") or []
return [k for k in [MPIEXEC,MPIEXEC_NUMPROC_FLAG,str(nproc)] \
......
# File for module specific CMake tests.
# find all required packages
FIND_PACKAGE (HDF5 REQUIRED COMPONENTS C)
if(HDF5_IS_PARALLEL)
message(STATUS "Parallel HDF5 library found")
add_definitions(-DHDF5_PARALLEL)
FIND_PACKAGE (HDF5 REQUIRED)
if(NOT HDF5_IS_PARALLEL)
message(SEND_ERROR "Parallel HDF5 must be installed!")
endif()
add_definitions(-DHDF5_PARALLEL)
FIND_PACKAGE (FFTW REQUIRED)
FIND_PACKAGE (SuperLU REQUIRED)
FIND_PACKAGE (MPI REQUIRED)
if(NOT dune-uggrid_FOUND)
MESSAGE(STATUS "Module dune-uggrid was not found. Falling back to external UG grid manager")
FIND_PACKAGE(UG REQUIRED)
endif()
FIND_PACKAGE (METIS)
FIND_PACKAGE (ParMETIS)
......
function(scrape_parameters SOURCE_DIR XML_FILE CSS OUTPUT RESULT_NAME)
message(STATUS "Running parameter scraper on sources ${SOURCE_DIR}")
if(DEPLOY_SPHINX_SOURCE_URL)
execute_process(COMMAND ${CMAKE_BINARY_DIR}/dune-env scrape_folder.py --source ${SOURCE_DIR}
execute_process(COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env scrape_folder.py --source ${SOURCE_DIR}
--xml ${XML_FILE} --out ${OUTPUT} --css ${CSS}
--source_url ${DEPLOY_SPHINX_SOURCE_URL}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
RESULT_VARIABLE RETURN_CODE)
else()
execute_process(COMMAND ${CMAKE_BINARY_DIR}/dune-env scrape_folder.py --source ${SOURCE_DIR}
execute_process(COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env scrape_folder.py --source ${SOURCE_DIR}
--xml ${XML_FILE} --out ${OUTPUT} --css ${CSS}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
RESULT_VARIABLE RETURN_CODE)
......@@ -18,11 +18,11 @@ endfunction()
file(COPY ${CMAKE_CURRENT_LIST_DIR} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/..)
scrape_parameters(${PROJECT_SOURCE_DIR}/dune/dorie-rfg ${CMAKE_CURRENT_SOURCE_DIR}/field-parameters.xml ${CMAKE_CURRENT_SOURCE_DIR}/parameters.css "parfield.ini;field-parameters.html;field-parameters.rst" FIELDPARSCRAPE_RETURN)
if (${FIELDPARSCRAPE_RETURN} GREATER 0)
if (NOT ${FIELDPARSCRAPE_RETURN} EQUAL 0)
message(FATAL_ERROR "Parameter scraper failed. DORiE can not be built.")
endif()
scrape_parameters(${PROJECT_SOURCE_DIR}/dune/dorie ${CMAKE_CURRENT_SOURCE_DIR}/parameters.xml ${CMAKE_CURRENT_SOURCE_DIR}/parameters.css "config.ini;parameters.html;parameters.rst" PARSCRAPE_RETURN)
if (${PARSCRAPE_RETURN} GREATER 0)
if (NOT ${PARSCRAPE_RETURN} EQUAL 0)
message(FATAL_ERROR "Parameter scraper failed. DORiE can not be built.")
endif()
......@@ -7,5 +7,5 @@ Module: dorie
Version: 0.9
Maintainer: dorieteam@iup.uni-heidelberg.de
#depending on
Depends: dune-python dune-pdelab dune-randomfield
Suggests: dune-testtools dune-uggrid
Depends: dune-pdelab dune-uggrid dune-randomfield
Suggests: dune-testtools
\ No newline at end of file
......@@ -46,8 +46,10 @@
#include <dune/pdelab/common/function.hh>
#include <dune/pdelab/constraints/p0.hh>
#include <dune/pdelab/constraints/p0ghost.hh>
#include <dune/pdelab/constraints/conforming.hh>
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/finiteelementmap/qkdg.hh>
#include <dune/pdelab/finiteelementmap/p0fem.hh>
#include <dune/pdelab/finiteelementmap/pkfem.hh>
#include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include <dune/pdelab/gridfunctionspace/localfunctionspace.hh>
......
......@@ -11,17 +11,18 @@ Simulation<Traits>::Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid>
verbose(inifile.get<int>("output.verbose"))
{
Dune::Timer timer;
const int verbose_low = helper.rank() == 0 ? verbose : 0;
// --- 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);
Dune::PDELab::constraints(*gfs,*cc,false);
// --- Parameter class ---
InterpolatorFactory intfac(inifile,helper,verbose);
InterpolatorFactory intfac(inifile,helper,verbose_low);
interp = intfac.create();
ParameterizationFactory pfac(inifile,helper,*interp,verbose);
ParameterizationFactory pfac(inifile,helper,*interp,verbose_low);
param = pfac.create();
// --- Operator Helper Classes ---
......@@ -63,20 +64,44 @@ void Simulation<Traits>::operator_setup ()
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);
// build parallel solvers
if (helper.size() > 1) {
// --- Solvers ---
lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
lscc = std::make_unique<LSCC>();
ls_par = std::make_unique<LSPar>(*igo,*cc,*lsgfs,*lscc,1000,0,true,true);
pdesolver_par = std::make_unique<PDESOLVERPar>(*igo,*ls_par);
pdesolver_par->setParameters(inifile.sub("NewtonParameters"));
pdesolver_par->setVerbosityLevel(verbose);
// --- Time Step Operators ---
Dune::PDELab::OneStepMethodResult osm_result;
if(osm_par){
osm_result = osm_par->result(); // cache old result
}
osm_par = std::make_unique<OSMPar>(alex2,*igo,*pdesolver_par);
osm_par->setResult(osm_result);
osm_par->setVerbosityLevel(verbose+1);
}
// build sequential solvers
else {
ls_seq = std::make_unique<LSSeq>(5000,0);
pdesolver_seq = std::make_unique<PDESOLVERSeq>(*igo,*ls_seq);
pdesolver_seq->setParameters(inifile.sub("NewtonParameters"));
pdesolver_seq->setVerbosityLevel(verbose);
// --- Time Step Operators ---
Dune::PDELab::OneStepMethodResult osm_result;
if(osm){
osm_result = osm->result(); // cache old result
// --- Time Step Operators ---
Dune::PDELab::OneStepMethodResult osm_result;
if(osm_seq){
osm_result = osm_seq->result(); // cache old result
}
osm_seq = std::make_unique<OSMSeq>(alex2,*igo,*pdesolver_seq);
osm_seq->setResult(osm_result);
osm_seq->setVerbosityLevel(verbose+1);
}
osm = std::make_unique<OSM>(alex2,*igo,*pdesolver);
osm->setResult(osm_result);
osm->setVerbosityLevel(verbose+1);
gfs->update();
......@@ -102,11 +127,22 @@ 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);
// solve in parallel
if (helper.size() > 1) {
pdesolver_par->setMaxIterations(controller->getIterations());
dwarn.push(false);
osm_par->apply(t,dt,*uold,*unew);
dwarn.pop();
}
// solve sequentially
else {
pdesolver_seq->setMaxIterations(controller->getIterations());
osm_seq->apply(t,dt,*uold,*unew);
}
*uold = *unew;
}
catch (Dune::PDELab::NewtonError &e){
......
......@@ -23,6 +23,14 @@ protected:
/// GFS Constraints Container
using CC = typename GFSHelper::CC;
/// LSGFS Helper
using LSGFSHelper = Dune::Dorie::LinearSolverGridFunctionSpaceHelper<GV,RF,Traits::GridGeometryType>;
/// Linear solver GFS
using LSGFS = typename LSGFSHelper::Type;
/// LSGFS Constraints Container
using LSCC = typename LSGFSHelper::CC;
// -- DORiE Class Definitions -- //
/// Parameter Interpolator
using Interpolator = Dune::Dorie::InterpolatorBase<Traits,Traits::dim>;
......@@ -56,14 +64,18 @@ protected:
using IGO = Dune::PDELab::OneStepGridOperator<GO0,GO1>;
/// Solution vector type
using U = typename IGO::Traits::Domain;
/// Linear solver type
using LS = Dune::PDELab::ISTLBackend_SEQ_SuperLU;
/// Non-linear solver type
using PDESOLVER = Dune::PDELab::Newton<IGO,LS,U>;
/// Linear solver types
using LSSeq = Dune::PDELab::ISTLBackend_SEQ_SuperLU;
using LSPar = Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<IGO,CC,LSGFS,LSCC,
Dune::PDELab::CG2DGProlongation,Dune::SeqSSOR,Dune::BiCGSTABSolver>;
/// Non-linear solver types
using PDESOLVERSeq = Dune::PDELab::Newton<IGO,LSSeq,U>;
using PDESOLVERPar = Dune::PDELab::Newton<IGO,LSPar,U>;
/// Time stepping scheme
using TimeStepScheme = Dune::PDELab::Alexander2Parameter<RF>;
/// Method computing the time step
using OSM = Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVER,U,U>;
/// Methods computing the time step
using OSMSeq = Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVERSeq,U,U>;
using OSMPar = Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVERPar,U,U>;
// -- Utility Class Definitions -- //
/// VTK Output writer base class
......@@ -82,6 +94,9 @@ protected:
std::unique_ptr<GFS> gfs;
std::unique_ptr<CC> cc;
std::unique_ptr<LSGFS> lsgfs;
std::unique_ptr<LSCC> lscc;
std::unique_ptr<Interpolator> interp;
std::unique_ptr<Parameters> param;
std::unique_ptr<FlowBoundary> fboundary;
......@@ -96,9 +111,12 @@ protected:
std::unique_ptr<GO0> go0;
std::unique_ptr<GO1> go1;
std::unique_ptr<IGO> igo;
std::unique_ptr<LS> ls;
std::unique_ptr<PDESOLVER> pdesolver;
std::unique_ptr<OSM> osm;
std::unique_ptr<LSSeq> ls_seq;
std::unique_ptr<LSPar> ls_par;
std::unique_ptr<PDESOLVERSeq> pdesolver_seq;
std::unique_ptr<PDESOLVERPar> pdesolver_par;
std::unique_ptr<OSMSeq> osm_seq;
std::unique_ptr<OSMPar> osm_par;
TimeStepScheme alex2;
std::unique_ptr<U> uold;
......
......@@ -31,6 +31,77 @@ R estimate_mbe_entries (const int dim, const Dune::GeometryType::BasicType geo)
return 1;
}
/// Provide types and construction of the GFS for the linear solver
template<typename GridView, typename RF, Dune::GeometryType::BasicType GeometryType>
struct LinearSolverGridFunctionSpaceHelper
{};
/// Provide types and construction of the GridFunctionSpace
template<typename GridView, typename RF>
struct LinearSolverGridFunctionSpaceHelper<GridView,RF,Dune::GeometryType::BasicType::simplex>
{
private:
static constexpr int dim = GridView::dimension;
using DF = typename GridView::Grid::ctype;
public:
/// Entity set of the GFS
using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
/// FiniteElementMap type of GFS
using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
/// Constraints type of the GFS
using CON = Dune::PDELab::P0ParallelConstraints;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed>>;
/// GFS Constraints Container type
using CC = typename Type::template ConstraintsContainer<RF>::Type;
/// create GFS from GridView
static Type create (const GridView& gv)
{
ES es(gv);
Dune::GeometryType geo;
geo.makeSimplex(dim);
auto fem = std::make_shared<FEM>(geo);
auto con = std::make_shared<CON>();
return Type(es,fem,con);
}
};
/// Provide types and construction of the GridFunctionSpace
template<typename GridView, typename RF>
struct LinearSolverGridFunctionSpaceHelper<GridView,RF,Dune::GeometryType::BasicType::cube>
{
private:
static constexpr int dim = GridView::dimension;
using DF = typename GridView::Grid::ctype;
public:
/// Entity set of the GFS
using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
/// FiniteElementMap type of GFS
using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
/// Constraints type of the GFS
using CON = Dune::PDELab::P0ParallelConstraints;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed> >;
/// GFS Constraints Container type
using CC = typename Type::template ConstraintsContainer<RF>::Type;
/// create GFS from GridView
static Type create (const GridView& gv)
{
ES es(gv);
Dune::GeometryType geo;
geo.makeCube(dim);
auto fem = std::make_shared<FEM>(geo);
auto con = std::make_shared<CON>();
return Type(es,fem,con);
}
};
/// Provide types and construction of the GridFunctionSpace. Has to be specialized.
template<typename GridView, typename RF, int order, Dune::GeometryType::BasicType GeometryType>
struct GridFunctionSpaceHelper
......@@ -46,21 +117,25 @@ private:
using DF = typename GridView::Grid::ctype;
public:
/// Entity set of the GFS
using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
/// FiniteElementMap type of GFS
using FEM = typename Dune::PDELab::PkLocalFiniteElementMap<GridView,DF,RF,order>;
using FEM = typename Dune::PDELab::PkLocalFiniteElementMap<ES,DF,RF,order>;
/// Constraints type of the GFS
using CON = Dune::PDELab::P0ParallelConstraints;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<GridView,FEM,
Dune::PDELab::P0ParallelConstraints,
Dune::PDELab::istl::VectorBackend<
Dune::PDELab::istl::Blocking::fixed> >;
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed>>;
/// GFS Constraints Container type
using CC = typename Type::template ConstraintsContainer<RF>::Type;
/// create GFS from GridView
static Type create (const GridView& gv)
{
auto fem = std::make_shared<FEM>(gv);
return Type(gv,fem);
ES es(gv);
auto fem = std::make_shared<FEM>(es);
auto con = std::make_shared<CON>();
return Type(es,fem,con);
}
};
......@@ -73,21 +148,25 @@ private:
using DF = typename GridView::Grid::ctype;
public:
/// Entity set of the GFS
using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
/// FiniteElementMap type of GFS
using FEM = typename Dune::PDELab::QkDGLocalFiniteElementMap<DF,RF,order,dim>;
/// Constraints type of the GFS
using CON = Dune::PDELab::P0ParallelConstraints;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<GridView,FEM,
Dune::PDELab::P0ParallelGhostConstraints,
Dune::PDELab::istl::VectorBackend<
Dune::PDELab::istl::Blocking::fixed> >;
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed> >;
/// GFS Constraints Container type
using CC = typename Type::template ConstraintsContainer<RF>::Type;
/// create GFS from GridView
static Type create (const GridView& gv)
{
ES es(gv);
auto fem = std::make_shared<FEM>();
return Type(gv,fem);
auto con = std::make_shared<CON>();
return Type(es,fem,con);
}
};
......
This diff is collapsed.
#ifndef DUNE_DORIE_ERROR_INDICATOR_HH
#define DUNE_DORIE_ERROR_INDICATOR_HH
#include<dune/geometry/referenceelements.hh>
#include<dune/pdelab/common/referenceelements.hh>
#include<dune/pdelab/common/quadraturerules.hh>
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
#include<dune/pdelab/localoperator/idefault.hh>
#include<dune/pdelab/localoperator/defaultimp.hh>
#include<dune/pdelab/finiteelement/localbasiscache.hh>
namespace Dune {
......@@ -85,23 +95,22 @@ namespace Dune {
// geometric factor of penalty
const RF h_F = std::min(ig.inside().geometry().volume(),ig.outside().geometry().volume())/ig.geometry().volume(); // Houston!
// select quadrature rule
Dune::GeometryType gtface = ig.geometry().type();
const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
// get element geometry
auto gtface = ig.geometry();
// container for sum of errors
RF sum(0.0);
int counter = 0;
// loop over quadrature points and integrate normal flux
for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it = rule.begin(); it != rule.end(); ++it)
const auto quadrature_rule = quadratureRule(gtface,intorder);
for (const auto& it : quadrature_rule)
{
const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it->position());
// retrieve unit normal vector
const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it.position());
// position of quadrature point in local coordinates of elements
const Domain local_s = ig.geometryInInside().global(it->position());
const Domain local_n = ig.geometryInOutside().global(it->position());
const Domain local_s = ig.geometryInInside().global(it.position());
const Domain local_n = ig.geometryInOutside().global(it.position());
// evaluate basis functions
std::vector<Scalar> phi_s(lfsu_s.size());
......@@ -142,15 +151,22 @@ namespace Dune {
for (unsigned int i = 0; i<lfsu_n.size(); i++)
gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
// offset for parameter queries
const RF eps = 1e-9;
// parameter queries: inside
const Domain xGlobal_s = ig.inside().geometry().global(local_s);
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);
// outside
const Domain xGlobal_n = ig.outside().geometry().global(local_n);
Domain xGlobal_n = ig.outside().geometry().global(local_n);
xGlobal_n = xGlobal_n.axpy(eps,n_F);
const RF head_n = param.headMillerToRef(u_n,xGlobal_n);
const RF satCond_n = param.condRefToMiller(param.K(xGlobal_n),xGlobal_n);
const RF saturation_n = param.saturation(head_n,xGlobal_n);
......@@ -182,14 +198,14 @@ namespace Dune {
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());
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
// add square of error
sum += jump*jump*factor*factor;
counter++;
}
// correct for the number of neighbors
sum *= 1. / counter;
sum *= 1. / quadrature_rule.size();
// what is this variable for?
DF h_T = 1;//std::max( diameter(ig.inside().geometry()),
......@@ -225,3 +241,5 @@ namespace Dune {
}
}
#endif // DUNE_DORIE_ERROR_INDICATOR_HH
\ No newline at end of file
......@@ -219,11 +219,8 @@ namespace Dune {
/// Opens H5 array for reading
void open_array(std::string groupName)
{
if(helper.rank()==0)
{
h5file.open(groupName);
array_open = true;
}
h5file.open(groupName);
array_open = true;
}
/// Reads parameter array for Parameter p from H5 file
......@@ -232,12 +229,9 @@ namespace Dune {
if (!array_open)
DUNE_THROW(Exception,"Array file must be opened before reading a parameter");
if(helper.rank()==0)
{
Array fieldArray;
h5file.read_dataset(fieldArray, fieldCells, p.name);
p.set_values(fieldArray);
}
Array fieldArray;
h5file.read_dataset(fieldArray, fieldCells, p.name);
p.set_values(fieldArray);
}
/// Closes H5 file
......@@ -246,11 +240,8 @@ namespace Dune {
if (!array_open)
DUNE_THROW(Exception,"Array file already closed");
if(helper.rank()==0)
{
h5file.close();
array_open = false;
}
h5file.close();
array_open = false;
}
};
......
......@@ -4,6 +4,54 @@
namespace Dune {
namespace Dorie {
/// Check if parallel execution is allowed
/** Currently, we only allow parallel execution on YaspGrid due to
* constraints imposed by the linear solver
*/
template<typename GridType>
void check_parallel_allowed (const Dune::MPIHelper& helper)
{
if(helper.size() > 1
&& !std::is_same<YaspGrid<GridType::dimension>,GridType>::value) {
DUNE_THROW(Dune::Exception,"DORiE does not support parallel execution with this grid configuration!");
}
}
/// Specialized constructor for parallel YaspGrid construction
/** \param extensions Extensions of the grid
* \param elements Number of elements in each direction
* \return Shared pointer to the grid
*/
template<typename GridType>
std::enable_if_t<
std::is_same<YaspGrid<GridType::dimension>,GridType>::value,
std::shared_ptr<YaspGrid<GridType::dimension>> >
grid_cube_construction (const Dune::FieldVector<double,GridType::dimension>& extensions,
const std::array<unsigned int,GridType::dimension>& elements)
{
std::array<int,GridType::dimension> elements_;
std::copy(elements.cbegin(),elements.cend(),elements_.begin());
return std::make_shared<YaspGrid<GridType::dimension>>(extensions,elements_);
}
/// Default rectangular grid constructor. Call StructuredGridFactory
/** \param extensions Extensions of the grid
* \param elements Number of elements in each direction
* \return Shared pointer to the grid
*/
template<typename GridType>
std::enable_if_t<
!std::is_same<YaspGrid<GridType::dimension>,GridType>::value,
std::shared_ptr<GridType> >
grid_cube_construction (const Dune::FieldVector<double,GridType::dimension>& extensions,
const std::array<unsigned int,GridType::dimension>& elements)
{
const Dune::FieldVector<double,GridType::dimension> origin(.0);
return Dune::StructuredGridFactory<GridType>::createCubeGrid(origin,extensions,elements);
}
/// Read simplex Gmsh file in 2D or 3D, create a grid, and return a shared pointer to the grid
/** Refine the grid globally and call default load balancing for multi-processing support.
* \tparam GridType Type of the grid. Must be supported by Dune::GridFactory
......@@ -13,13 +61,14 @@ namespace Dorie {
template<class GridType>
std::shared_ptr<GridType> build_grid_gmsh (const Dune::ParameterTree& inifile, const Dune::MPIHelper& helper)
{
check_parallel_allowed<GridType>(helper);
Dune::Timer timer;
ios_base_all_saver restore(std::cout);
enum { dim = GridType::dimension };
const int verbose = inifile.get<int>("output.verbose");
const int level = inifile.get<int>("grid.initialLevel");
const std::string meshfilename = inifile.get<std::string>("grid.gridFile");
typedef std::vector<int> GmshIndexMap;
if(verbose>0 && helper.rank()==0){
std::cout << "GRID SETUP:" << std::endl;
......@@ -29,14 +78,17 @@ std::shared_ptr<GridType> build_grid_gmsh (const Dune::ParameterTree& inifile, c
}
auto grid = std::make_shared<GridType>();
GmshIndexMap boundary_index_map;
GmshIndexMap element_index_map;
Dune::GridFactory<GridType> factory(grid.get());
Dune::GmshReader<GridType>::read(factory,meshfilename,boundary_index_map,element_index_map,verbose>2?true:false,false);
factory.createGrid();
grid->globalRefine(level);
if(helper.rank() == 0){
std::vector<int> boundary_index_map;
std::vector<int> element_index_map;
Dune::GmshReader<GridType>::read(factory,meshfilename,boundary_index_map,element_index_map,verbose>2?true:false,