Commit cdb7a395 authored by Santiago Ospina's avatar Santiago Ospina

Merge remote-tracking branch...

Merge remote-tracking branch 'origin/98-implement-a-simulation-object-for-transport-transportsimulation' into 101-couple-the-transportsimulation-with-richardssimulation
parents 911c397d e5c13de6
......@@ -165,4 +165,4 @@ create_default_config(
"config.ini;parameters.html;parameters.rst"
${PROJECT_SOURCE_DIR}/dune/dorie
${CMAKE_CURRENT_SOURCE_DIR}/parameters.css
)
)
\ No newline at end of file
add_subdirectory("groups")
# shortcut for creating the Doxyfile.in and Doxyfile
add_doxygen_target()
......@@ -3,6 +3,7 @@
# Location of source files
INPUT += @top_srcdir@/dune/dorie/
INPUT += @top_srcdir@/doc/doxygen/groups
# extract private class members
EXTRACT_PRIVATE = YES
......@@ -18,5 +19,4 @@ USE_MATHJAX = YES
# subdirectory from a directory tree whose root is specified with the INPUT tag.
# exclude versions of the operators that are not used
# EXCLUDE += @top_srcdir@/src/operator_DG_old.hh \
# @top_srcdir@/src/operator_DG_no_upwind.hh
# EXCLUDE += @top_srcdir@/...
/**
@defgroup Common Common
@{
@todo document models
@}
**/
\ No newline at end of file
/** @defgroup LocalOperators Local operators
@{
Local operators are in many senses the heart of dorie; in them, we transform
finite element basis into a vector of residuals and a mass matrix taking into
account the specific equation and method we want to solve. In particular,
some care has been taken regarding the different frames of reference when
working with coordinates and gradients. It can get particularly messy
particularly when working with intersections because one has to express the
same position/gradient with respect to the inside, outside, and face entities,
and sometimes with respect to the grid (e.g. global coordinates). Therefore,
in dorie we try to use the following suffixes convention:
| Convention | Meaning |
|------------|---------|
| `<description>` + `_i` | `<description>` with respect to the **inside** entity. |
| `<description>` + `_o` | `<description>` with respect to the **outside** entity. |
| `<description>` + `_f` | `<description>` with respect to the **face** entity (or intersection). |
| `<description>` + `_g` | `<description>` with respect to the **global** grid. |
| `<description>` + `_n` | `<description>` in **normal** direction. |
Notice that when working with gradients, the suffix for normal direction can
be combined with the first three suffixes (`_n_i`, `_n_o`, and `_n_f`). When it's
not specified, is understood that it is with respect to the face
(`_n` \f$\equiv\f$. `_n_f`).
@todo Update Dune::Dorie::Operator::RichardsDGSpatialOperator to the convention.
@}
**/
/**
@defgroup Models Models
@{
@todo document models
@}
@defgroup RichardsModel Richards
@{
@ingroup Models
@todo document richards model
@}
@defgroup TransportModel Transport
@{
@ingroup Models
@todo document transport model
@}
**/
\ No newline at end of file
......@@ -4,11 +4,12 @@ if(dune-testtools_FOUND)
add_subdirectory("test")
endif()
add_executable("transport" transport.cc)
dune_target_link_libraries(transport dorie-impl ${DUNE_LIBS})
add_executable("richards" richards.cc)
dune_target_link_libraries(richards dorie-impl ${DUNE_LIBS})
dune_target_link_libraries(richards richards-impl ${DUNE_LIBS})
add_executable("transport" transport.cc)
dune_target_link_libraries(transport richards-impl ${DUNE_LIBS})
dune_target_link_libraries(transport transport-impl ${DUNE_LIBS})
add_custom_target("dorie" DEPENDS richards transport)
......
#install headers
install(FILES
boundary_condition.hh
cfl_condition.hh
grid_creator.hh
grid_function_writer.hh
grid_mapper.hh
h5tools.hh
interpolator.hh
setup_inifile.hh
......
......@@ -36,6 +36,27 @@ namespace Setup
return richards_ini;
}
/*------------------------------------------------------------------------*//**
* @brief Set up the configuration keys for Transport simulations
*
* @param ini The inifile
*/
Dune::ParameterTree prep_ini_for_transport(const Dune::ParameterTree& ini)
{
Dune::ParameterTree richards_ini(ini.sub("transport"));
// move grid extensions into richards category (needed for reading parameters)
richards_ini["grid.extensions"] = ini["grid.extensions"];
const std::string adapt_policy_str = ini.get<std::string>("adaptivity.policy");
if (adapt_policy_str != "none") {
DUNE_THROW(NotImplemented, "Not known adaptativity policy ("
<< ini.get<std::string>("adaptivity.policy") << ")!");
}
return richards_ini;
}
};
}
......
add_subdirectory("richards")
\ No newline at end of file
add_subdirectory("richards")
add_subdirectory("transport")
\ No newline at end of file
......@@ -15,7 +15,6 @@
#include <dune/pdelab/localoperator/defaultimp.hh>
#include <dune/pdelab/finiteelement/localbasiscache.hh>
#include <dune/dorie/model/richards/local_operator.hh>
namespace Dune {
......
add_library(dorie-impl STATIC
add_library(richards-impl STATIC
sim_yasp_2_1.cc
sim_yasp_2_2.cc
sim_ug_2_1.cc
......
......@@ -79,7 +79,7 @@ inline auto read_experimental_operator_settings(
RichardsDGUpwinding::Type upwinding;
RichardsDGWeights::Type weights;
const auto method_str = inifile.get<std::string>("numerics.experimental.method");
const auto method_str = inifile.get<std::string>("dg.experimental.method");
if (method_str == "SIPG")
method = RichardsDGMethod::SIPG;
else if (method_str == "NIPG")
......@@ -92,7 +92,7 @@ inline auto read_experimental_operator_settings(
DUNE_THROW(Dune::IOError, "Unknown DG method '" + method_str + "'!");
const auto upwinding_str
= inifile.get<std::string>("numerics.experimental.upwinding");
= inifile.get<std::string>("dg.experimental.upwinding");
if (upwinding_str == "none")
upwinding = RichardsDGUpwinding::none;
else if (upwinding_str == "semiUpwind")
......@@ -102,7 +102,7 @@ inline auto read_experimental_operator_settings(
else
DUNE_THROW(Dune::IOError, "Unknown upwinding '" + upwinding_str + "'!");
const auto weights_str = inifile.get<bool>("numerics.experimental.weights");
const auto weights_str = inifile.get<bool>("dg.experimental.weights");
if (weights_str)
weights = RichardsDGWeights::weightsOn;
else
......@@ -122,6 +122,8 @@ inline auto read_experimental_operator_settings(
* to allow for different orders of the shape functions, which should be
* enough to support p-adaptivity. (Another likely candidate would be
* differing geometry types, i.e. hybrid meshes.)
*
* @ingroup LocalOperators
*/
template<typename Traits, typename Parameter, typename Boundary, typename SourceTerm, typename FEM, bool adjoint>
class RichardsDGSpatialOperator
......@@ -966,6 +968,8 @@ private:
/**
* @brief DG local operator for the temporal derivative of the Richards equation
*
* @ingroup LocalOperators
*/
template<typename Traits, typename Parameter, typename FEM, bool adjoint>
class RichardsDGTemporalOperator
......@@ -1071,4 +1075,4 @@ template<typename Traits, typename Parameter, typename FEM, bool adjoint>
} // namespace Dune
#endif // DUNE_DORIE_RICHARDS_OPERATOR_HH
#endif // DUNE_DORIE_RICHARDS_OPERATOR_HH
\ No newline at end of file
......@@ -34,7 +34,6 @@
#include <dune/dorie/model/richards/flow_source.hh>
#include <dune/dorie/model/richards/local_operator.hh>
namespace Dune{
namespace Dorie{
......@@ -43,6 +42,8 @@ namespace Dorie{
* class extends BaseTraits to be used for the DG implementation of
* the Richards simulation
*
* @ingroup RichardsModel
*
* @tparam BaseTraits Traits defining domain and range field properties of
* the simulation.
* @tparam order Order of the polynomial degree used for the basis
......@@ -140,6 +141,8 @@ struct RichardsSimulationTraits : public BaseTraits
* for the richards equation. It can perform time-steps and print
* the solution
*
* @ingroup RichardsModel
*
* @tparam RichardsSimulationTraits Traits containing the type definitions
* for Richards Simulation
*/
......
add_subdirectory("impl")
#install headers
install(FILES
solute_boundary.hh
solute_initial.hh
solute_source.hh
local_operator_FV.hh
transport.hh
richards_coupling.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/dorie/model/transport)
\ No newline at end of file
add_library(transport-impl STATIC
sim_yasp_2_0.cc
# sim_yasp_2_1.cc
# sim_yasp_2_2.cc
sim_ug_2_0.cc
# sim_ug_2_1.cc
# sim_yasp_2_3.cc
sim_yasp_3_0.cc
# sim_yasp_3_1.cc
# sim_ug_2_2.cc
# sim_yasp_3_2.cc
# sim_yasp_3_3.cc
# sim_ug_2_3.cc
sim_ug_3_0.cc
# sim_ug_3_1.cc
# sim_ug_3_2.cc
# sim_ug_3_3.cc
)
\ No newline at end of file
#ifndef DUNE_DORIE_TRANSPORT_IMPL_HH
#define DUNE_DORIE_TRANSPORT_IMPL_HH
#include <dune/grid/yaspgrid.hh>
#include <dune/geometry/type.hh>
namespace Dune{
namespace Dorie{
using Geo = Dune::GeometryType::BasicType;
}
}
#endif // DUNE_DORIE_TRANSPORT_IMPL_HH
\ No newline at end of file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/model/transport/impl/impl.hh>
#include <dune/dorie/model/transport/transport.cc>
#include <dune/dorie/model/transport/richards_coupling.cc>
namespace Dune{
namespace Dorie{
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1,0>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2,0>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3,0>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1,0>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2,0>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/model/transport/impl/impl.hh>
#include <dune/dorie/model/transport/transport.cc>
#include <dune/dorie/model/transport/richards_coupling.cc>
namespace Dune{
namespace Dorie{
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,1,0>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,2,0>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,3,0>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,1,0>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,2,0>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,3,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/model/transport/impl/impl.hh>
#include <dune/dorie/model/transport/transport.cc>
#include <dune/dorie/model/transport/richards_coupling.cc>
namespace Dune{
namespace Dorie{
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,0>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,0>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/model/transport/impl/impl.hh>
#include <dune/dorie/model/transport/transport.cc>
#include <dune/dorie/model/transport/richards_coupling.cc>
namespace Dune{
namespace Dorie{
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,1>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,1>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,1>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/model/transport/impl/impl.hh>
#include <dune/dorie/model/transport/transport.cc>
#include <dune/dorie/model/transport/richards_coupling.cc>
namespace Dune{
namespace Dorie{
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,2>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,2>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,2>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/model/transport/impl/impl.hh>
#include <dune/dorie/model/transport/transport.cc>
#include <dune/dorie/model/transport/richards_coupling.cc>
namespace Dune{
namespace Dorie{
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,1,3>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2,3>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/model/transport/impl/impl.hh>
#include <dune/dorie/model/transport/transport.cc>
#include <dune/dorie/model/transport/richards_coupling.cc>
namespace Dune{
namespace Dorie{
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1,0>>;
template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2,0>>;
// template class RichardsTransportCouplingSimulation<RichardsTransportCouplingSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3,0>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -12,7 +12,7 @@
#include <dune/pdelab/localoperator/idefault.hh>
#include <dune/pdelab/localoperator/defaultimp.hh>
#include "util_typedefs.hh"
#include <dune/dorie/common/typedefs.hh>
namespace Dune {
namespace Dorie {
......@@ -21,6 +21,7 @@ namespace Operator {
/*-------------------------------------------------------------------------*//**
* @brief Spatial local operator for the transport equation in unsaturated
* media in a finite volume scheme.
*
* @f{eqnarray*}{
* \partial_t[\theta C_w] +
* \nabla\cdot [\textbf{j}_w C_w] -
......@@ -33,37 +34,40 @@ namespace Operator {
* \Gamma_D
* @f}
*
*
* @todo Use diffusion coefficient grid function.
* @todo Implement outflow boundary condition.
*
* @ingroup LocalOperators
*
* @tparam Boundary Type of the class providing boundary
* conditions
* @tparam GridFunctionWaterFlux Type of a grid function which provides
* the water flux
* @tparam GridFunctionSaturation Type of a grid function which provides
* the saturation
* @tparam Boundary Type of the class providing boundary
* conditions
* @tparam GridFunctionWaterFlux Type of a grid function which provides
* the water flux
* @tparam GridFunctionWaterContent Type of a grid function which provides
* the water content
*/
template<class Boundary, class GridFunctionWaterFlux, class GridFunctionSaturation>
template<class Boundary, class GridFunctionWaterFlux, class GridFunctionWaterContent>
class TransportFVSpatialOperator
: public Dune::PDELab::NumericalJacobianVolume<
TransportFVSpatialOperator<
Boundary,
GridFunctionWaterFlux,
GridFunctionSaturation
GridFunctionWaterContent
>
>
, public Dune::PDELab::NumericalJacobianSkeleton<
TransportFVSpatialOperator<
Boundary,
GridFunctionWaterFlux,
GridFunctionSaturation
GridFunctionWaterContent
>
>
, public Dune::PDELab::NumericalJacobianBoundary<
TransportFVSpatialOperator<
Boundary,
GridFunctionWaterFlux,
GridFunctionSaturation
GridFunctionWaterContent
>
>
, public Dune::PDELab::FullSkeletonPattern
......@@ -85,37 +89,31 @@ public:
private:
static_assert(std::is_same<
typename GridFunctionWaterFlux::Traits::GridViewType,
typename GridFunctionSaturation::Traits::GridViewType>::value,
typename GridFunctionWaterContent::Traits::GridViewType>::value,
"TransportFVSpatialOperator: GridFunctionWaterFlux and"
"GridFunctionSaturation has to use the same grid view.");
// Define domain field
using DF = typename GridFunctionSaturation::Traits::DomainFieldType;
"GridFunctionWaterContent has to use the same grid view.");
using WaterFlux =typename GridFunctionWaterFlux::Traits::RangeType;
using Saturation =typename GridFunctionSaturation::Traits::RangeType;
// Extract world dimension
enum {dim = GridFunctionSaturation::Traits::dimDomain};
using WaterContent =typename GridFunctionWaterContent::Traits::RangeType;
public:
/**
* @brief Constructor of TransportFVSpatialOperator
*
* @param boundary The boundary terms
* @param[in] gf_flux The grid function of water flux
* @param[in] gf_sat The grid function of water saturation
* @param[in] diff_coeff The diffusion coefficient
* @param boundary The boundary terms
* @param[in] gf_water_flux The grid function of water flux
* @param[in] gf_water_content The grid function of water content
* @param[in] diff_coeff The diffusion coefficient
*/
TransportFVSpatialOperator(
Boundary& boundary,
std::shared_ptr<GridFunctionWaterFlux> gf_flux,
std::shared_ptr<GridFunctionSaturation> gf_sat,
std::shared_ptr<GridFunctionWaterFlux> gf_water_flux,
std::shared_ptr<GridFunctionWaterContent> gf_water_content,
double diff_coeff
) : _boundary(boundary)
, _gf_flux(gf_flux)
, _gf_sat(gf_sat)
, _gf_water_flux(gf_water_flux)
, _gf_water_content(gf_water_content)
, _diff_coeff(diff_coeff)
{}
......@@ -166,43 +164,66 @@ public:
const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o,
R& r_i, R& r_o) const
{
// Define range field
using RF = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
// Get entity entities from both sides of the intersection
// Get local basis traits from local function space
using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits;
// Get range for trial space
using RangeU = typename LocalBasisTraitsU::RangeType;
// Assert that polynomial degree is always 0
assert(lfsu_i.finiteElement().localBasis().order()==0);
assert(lfsu_o.finiteElement().localBasis().order()==0);
assert(lfsv_i.finiteElement().localBasis().order()==0);
assert(lfsv_o.finiteElement().localBasis().order()==0);
// Define entities for each frame of reference
const auto& entity_f = ig.intersection();
const auto& entity_i = ig.inside();
const auto& entity_o = ig.outside();
// Get geometries
auto geo = ig.geometry();
auto geo_f = entity_f.geometry();
auto geo_i = entity_i.geometry();
auto geo_o = entity_o.geometry();
// Get geometry of intersection in local coord. of neighbor entities
auto geo_in_i = ig.geometryInInside();
auto geo_in_i = entity_f.geometryInInside();
auto geo_in_o = entity_f.geometryInOutside();
// Face volume for integration
auto face_volume = geo.volume();
// Get volume of entities
const auto volume_f = geo_f.volume();
const auto volume_i = geo_i.volume();
const auto volume_o = geo_o.volume();
// Entity centers in references elements
const auto& ref_el_i = referenceElement(geo_i);
const auto& ref_el_o = referenceElement(geo_o);
const auto& center_position_i = ref_el_i.position(0,0);
const auto& center_position_o = ref_el_o.position(0,0);
// Evaluate water content
WaterContent water_content_i,water_content_o;
_gf_water_content->evaluate(entity_i,center_position_i,water_content_i);
_gf_water_content->evaluate(entity_o,center_position_o,water_content_o);
// Evaluate diffusion coeff. from both sides and take harmonic average
auto diff_coeff_i = _diff_coeff; // FIXME: evaluate from grid function
auto diff_coeff_o = _diff_coeff; // FIXME: evaluate from grid function
auto normal_face = ig.centerUnitOuterNormal();
auto diff_coeff_avg = 2.0/( 1.0/(diff_coeff_i+ 1E-30)
+ 1.0/(diff_coeff_o+ 1E-30));
auto diff_coeff_i = _diff_coeff * water_content_i; // FIXME: evaluate from grid function
auto diff_coeff_o = _diff_coeff * water_content_o; // FIXME: evaluate from grid function
auto normal_f = entity_f.centerUnitOuterNormal();
auto diff_coeff_f = 2.0/( 1.0/(diff_coeff_i + 1E-30)
+ 1.0/(diff_coeff_o + 1E-30) );
// Get center position of the face w.r.t inside entity
auto p_center_face_i = geo_in_i.center();
auto face_position_i = geo_in_i.center();
// Evaluate convective term (i.e water flux)
// Evaluate velocity field assume H(div) velocity field => may choose any side
WaterFlux water_flux_i;
_gf_flux->evaluate(entity_i, p_center_face_i, water_flux_i);
auto water_flux_n = water_flux_i*normal_face;
_gf_water_flux->evaluate(entity_i, face_position_i, water_flux_i);
auto water_flux_n = water_flux_i*normal_f;
// Inside/outside solute value
RF u_i = x_i(lfsu_i,0);
RF u_o = x_o(lfsu_o,0);
RangeU u_i = x_i(lfsu_i,0);
RangeU u_o = x_o(lfsu_o,0);