Commit 95d86194 authored by Santiago Ospina's avatar Santiago Ospina

Merge branch 'master' into 98-implement-a-simulation-object-for-transport-transportsimulation

# Conflicts:
#	doc/default_files/CMakeLists.txt
#	dune/dorie/model/richards/adaptivity/operator_error_indicator.hh
#	dune/dorie/model/richards/local_operator.hh
#	dune/dorie/model/richards/richards.hh
#	dune/dorie/solver/operator_DG.hh
#	dune/dorie/solver/richards_operator_DG.hh
Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent b22b4e61
/**
@defgroup Common Common
@{
@todo document models
@}
**/
\ No newline at end of file
......@@ -23,5 +23,6 @@
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
@todo document models
@}
@defgroup RichardsModel Richards
@ingroup Models
@{
@todo document richards model
@ingroup Models
@todo document richards model
@}
@defgroup TransportModel Transport
@ingroup Models
@{
@todo document transport model
@ingroup Models
@todo document transport model
@}
**/
\ No newline at end of file
......@@ -40,34 +40,34 @@ namespace Operator {
*
* @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
......@@ -89,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)
{}
......@@ -170,74 +164,87 @@ 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);
// Upwinding
const auto& u_upwind = (water_flux_n>=0) ? u_i : u_o;
// Entity centers in references elements
const auto& ref_el_i = referenceElement(geo_i);
const auto& ref_el_o = referenceElement(geo_o);
const auto& p_center_i = ref_el_i.position(0,0);
const auto& p_center_o = ref_el_o.position(0,0);
// Entity centers in global coordinates
auto p_center_i_global = geo_i.global(p_center_i);
auto p_center_o_global = geo_o.global(p_center_o);
auto center_position_i_g = geo_i.global(center_position_i);
auto center_position_o_g = geo_o.global(center_position_o);
// Distance between the two entity centers
p_center_i_global -= p_center_o_global;
auto distance = p_center_i_global.two_norm();
center_position_i_g -= center_position_o_g;
auto distance = center_position_i_g.two_norm();
// Finite difference of u between the two entities
RF dudn = (u_o-u_i)/distance;
// Evaluate saturation
typename GridFunctionSaturation::Traits::RangeType sat_i;
_gf_sat->evaluate(entity_i,p_center_i,sat_i);
RangeU dudn = (u_o-u_i)/distance;
// Solute flux in normal direction w.r.t the intersection
auto normal_flux = u_upwind*water_flux_n - diff_coeff_avg*dudn*sat_i;
auto soulte_flux_n = u_upwind*water_flux_n - diff_coeff_f*dudn;
// Symmetric contribution to residual on inside element
r_i.accumulate(lfsv_i, 0, normal_flux*face_volume );
r_o.accumulate(lfsv_o, 0, -normal_flux*face_volume );
r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f );
r_o.accumulate(lfsv_o, 0, -soulte_flux_n*volume_f );
}
......@@ -261,66 +268,72 @@ public:
void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i,
const LFSV& lfsv_i, R& r_i) const
{
// 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(lfsv_i.finiteElement().localBasis().order()==0);
// Get entities
const auto& entitiy_face = ig.intersection();
const auto& entity_f = ig.intersection();
const auto& entity_i = ig.inside();
// Get geometries
auto geo = entitiy_face.geometry();
const auto& ref_el = referenceElement(geo);
const auto& p_center_face = ref_el.position(0,0);
auto geo_f = entity_f.geometry();
const auto& ref_el_f = referenceElement(geo_f);
const auto& center_position_f = ref_el_f.position(0,0);
auto bc = _boundary.bc(entitiy_face,p_center_face,_time);
auto bc = _boundary.bc(entity_f,center_position_f,_time);
if (BoundaryCondition::isDirichlet(bc))
{
// Define range field
using RF = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
auto geo_i = entity_i.geometry();
auto p_center_i_global = geo_i.center();
auto center_position_i_g = geo_i.center();
// Face center in global coordinates
auto p_center_face_global = geo.center();
auto center_position_f_g = geo_f.center();
// Compute distance of these two points
p_center_i_global -= p_center_face_global;
auto distance = p_center_i_global.two_norm();
center_position_i_g -= center_position_f_g;
auto distance = center_position_i_g.two_norm();
// Evaluate Dirichlet condition
auto g = _boundary.g(entitiy_face,p_center_face,_time);
auto g = _boundary.g(entity_f,center_position_f,_time);
// Face volume for integration
auto face_volume = geo.volume();
auto volume_f = geo_f.volume();
// Get geometry of intersection in local coord. of neighbor entity
auto geo_in_i = ig.geometryInInside();
auto p_center_face_i = geo_in_i.center();
auto face_position_i = geo_in_i.center();
// Evaluate convective term (i.e water flux)
WaterFlux water_flux_i;
_gf_flux->evaluate(entity_i, p_center_face_i, water_flux_i);
auto normal_face = ig.unitOuterNormal(p_center_face);
auto water_flux_n = water_flux_i*normal_face;
_gf_water_flux->evaluate(entity_i, face_position_i, water_flux_i);
auto normal_f = ig.unitOuterNormal(center_position_f);
auto water_flux_n = water_flux_i*normal_f;
// Evaluate saturation
typename GridFunctionSaturation::Traits::RangeType sat_i;
_gf_sat->evaluate(entity_i, p_center_face_i, sat_i);
// Evaluate water content
WaterContent water_content_i;
_gf_water_content->evaluate(entity_i, face_position_i, water_content_i);
// Evaluate diffusion coefficient from inside
auto diff_coeff_i = _diff_coeff; // FIXME: evaluate from grid function
// Inside unknown value
RF u = x_i(lfsu_i,0);
RangeU u = x_i(lfsu_i,0);
// Only allow dirichlet values for influx cases
if (water_flux_n<=0){
// Solute flux in normal direction w.r.t the intersection
auto normal_flux = (g*water_flux_n)-(diff_coeff_i*sat_i*(g-u)/distance);
auto soulte_flux_n = (g*water_flux_n)-(diff_coeff_i*water_content_i*(g-u)/distance);
// Contribution to residual from Dirichlet boundary
r_i.accumulate(lfsv_i, 0, normal_flux*face_volume);
r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f);
}
......@@ -329,25 +342,25 @@ public:
{
// Get geometry of intersection in local coord. of neighbor entity
auto geo_in_i = ig.geometryInInside();
auto p_center_face_i = geo_in_i.center();
auto face_position_i = geo_in_i.center();
// Evaluate saturation
typename GridFunctionSaturation::Traits::RangeType sat_i;
_gf_sat->evaluate(entity_i, p_center_face_i, sat_i);
// Evaluate water content
typename GridFunctionWaterContent::Traits::RangeType water_content_i;
_gf_water_content->evaluate(entity_i, face_position_i, water_content_i);
// Face volume for integration
auto face_volume = geo.volume();
auto volume_f = geo_f.volume();
// Evaluate diffusion coefficient from inside
auto diff_coeff_i = _diff_coeff;
// Contribution to residual from Neumann boundary
auto j = _boundary.j(entitiy_face,p_center_face,_time);
auto j = _boundary.j(entity_f,center_position_f,_time);
// Solute flux in normal direction w.r.t the intersection
auto normal_flux = -diff_coeff_i*j*sat_i;
auto soulte_flux_n = -diff_coeff_i*j*water_content_i;
r_i.accumulate(lfsv_i, 0, normal_flux*face_volume );
r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f );
}
}
......@@ -397,8 +410,8 @@ public:
private:
Boundary& _boundary;
const std::shared_ptr<GridFunctionWaterFlux> _gf_flux;
const std::shared_ptr<GridFunctionSaturation> _gf_sat;
const std::shared_ptr<GridFunctionWaterFlux> _gf_water_flux;
const std::shared_ptr<GridFunctionWaterContent> _gf_water_content;
double _time;
double _diff_coeff;
};
......@@ -407,17 +420,17 @@ private:
* @brief Temporal local operator for the transport equation in unsaturated
* media in a finite volume scheme.
*
* @ingroup LocalOperators
*
* @tparam GridFunctionSaturation Type of a grid function which provides
* the saturation
* @ingroup LocalOperators
*
* @tparam GridFunctionWaterContent Type of a grid function which provides
* the water content
*/
template<class GridFunctionSaturation>
template<class GridFunctionWaterContent>
class TransportFVTemporalOperator
: public Dune::PDELab::NumericalJacobianVolume<
TransportFVTemporalOperator<GridFunctionSaturation>>
TransportFVTemporalOperator<GridFunctionWaterContent>>
, public Dune::PDELab::NumericalJacobianApplyVolume<
TransportFVTemporalOperator<GridFunctionSaturation>>
TransportFVTemporalOperator<GridFunctionWaterContent>>
, public Dune::PDELab::FullSkeletonPattern
, public Dune::PDELab::FullVolumePattern
, public Dune::PDELab::LocalOperatorDefaultFlags
......@@ -434,8 +447,8 @@ public:
public:
TransportFVTemporalOperator(
std::shared_ptr<GridFunctionSaturation> gfSaturation
) : _gfSaturation(gfSaturation)
std::shared_ptr<GridFunctionWaterContent> gf_water_content
) : _gf_water_content(gf_water_content)
{}
/*---------------------------------------------------------------------*//**
......@@ -457,25 +470,27 @@ public:
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
const LFSV& lfsv, R& r) const
{
// Define range field
using RF = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
// 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;
RF u = x(lfsu,0);
RangeU u = x(lfsu,0);
// entity geometry
auto geo = eg.geometry();
// inside cell center
const auto& ref_el = referenceElement(geo);
const auto& p_center = ref_el.position(0,0);
const auto& center_position = ref_el.position(0,0);
typename GridFunctionSaturation::Traits::RangeType saturation;
typename GridFunctionWaterContent::Traits::RangeType water_content;
_gfSaturation->evaluate(eg.entity(),p_center,saturation);
_gf_water_content->evaluate(eg.entity(),center_position,water_content);
// update residual
r.accumulate(lfsv ,0 , saturation*u*geo.volume());
r.accumulate(lfsv ,0 , water_content*u*geo.volume());
}
/*---------------------------------------------------------------------*//**
......@@ -492,7 +507,7 @@ public:
}
private:
std::shared_ptr<GridFunctionSaturation> _gfSaturation;
std::shared_ptr<GridFunctionWaterContent> _gf_water_content;
double _time;
};
......
......@@ -9,7 +9,7 @@ except NameError:
pass
# paths set by cmake
DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Debug/dorie"
DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Release/dorie"
MPIEXEC = "/usr/bin/mpiexec"
MPIEXEC_NUMPROC_FLAG = "-n"
MPIEXEC_PREFLAG = ""
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment