Commit 2a361e79 authored by Santiago Ospina's avatar Santiago Ospina
Browse files

Added constant initial condition, added todo-list to doxygen, and using upwinding always


Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 347410e8
...@@ -106,7 +106,11 @@ adding an empty line, make text **bold** or ``monospaced``. ...@@ -106,7 +106,11 @@ adding an empty line, make text **bold** or ``monospaced``.
</category> </category>
<category name="transport.initial"> <category name="transport.initial">
<!-- TODO! --> <parameter name="value">
<definition> Initial constant value over the whole domain. </definition>
<values> float </values>
<suggestion> 0 </suggestion>
</parameter>
</category> </category>
<category name="transport.time"> <category name="transport.time">
......
...@@ -27,11 +27,12 @@ TransportSimulation<Traits>::TransportSimulation( ...@@ -27,11 +27,12 @@ TransportSimulation<Traits>::TransportSimulation(
// --- Operator Helper Classes --- // --- Operator Helper Classes ---
sboundary = std::make_unique<SoluteBoundary>(inifile); sboundary = std::make_unique<SoluteBoundary>(inifile);
// ssource = std::make_unique<SoluteSource>(inifile,*param);
sinitial = std::make_unique<SoluteInitial>(gv,0.); sinitial = std::make_unique<SoluteInitial>(gv,0.);
// --- Solution Vectors and Initial Condition --- // --- Solution Vectors and Initial Condition ---
u = std::make_shared<U>(*gfs,.0); const RF ini_value = inifile.get<RF>("initial.value",0.);
u = std::make_shared<U>(*gfs,ini_value);
Dune::PDELab::interpolate(*sinitial,*gfs,*u); Dune::PDELab::interpolate(*sinitial,*gfs,*u);
controller = std::make_unique<CalculationController>(inifile,*sboundary,helper); controller = std::make_unique<CalculationController>(inifile,*sboundary,helper);
......
...@@ -34,10 +34,10 @@ namespace Dorie{ ...@@ -34,10 +34,10 @@ namespace Dorie{
* class extends BaseTraits to be used for the DG implementation of * class extends BaseTraits to be used for the DG implementation of
* the Transport simulation * the Transport simulation
* *
* @tparam BaseTraits Traits defining domain and range field * @tparam BaseTraits Traits defining domain and range field
* properties of the simulation. * properties of the simulation.
* @tparam GFSaturationType Grid Function type of the water content * @tparam GFWaterFluxType Grid Function type of the water flux
* @tparam GFWaterFluxType Grid Function type of the water flux * @tparam GFSaturationType Grid Function type of the saturation
*/ */
template<class BaseTraits, class GFWaterFluxType, class GFSaturationType> template<class BaseTraits, class GFWaterFluxType, class GFSaturationType>
struct TransportSimulationTraits : public BaseTraits struct TransportSimulationTraits : public BaseTraits
...@@ -119,6 +119,11 @@ struct TransportSimulationTraits : public BaseTraits ...@@ -119,6 +119,11 @@ struct TransportSimulationTraits : public BaseTraits
* for the transport equation. It can perform time-steps and print * for the transport equation. It can perform time-steps and print
* the solution. * the solution.
* *
* @todo Implement source term.
* @todo Implement more complex initial conditions.
* @todo Implement adaptivity.
* @todo Implement DG local operator.
*
* @tparam TransportSimulationTraits Traits containing the type definitions * @tparam TransportSimulationTraits Traits containing the type definitions
* which this class will use. * which this class will use.
*/ */
......
...@@ -2,14 +2,8 @@ ...@@ -2,14 +2,8 @@
#ifndef DUNE_DORIE_TRANSPORT_OPERATOR_HH #ifndef DUNE_DORIE_TRANSPORT_OPERATOR_HH
#define DUNE_DORIE_TRANSPORT_OPERATOR_HH #define DUNE_DORIE_TRANSPORT_OPERATOR_HH
#include <string>
#include <tuple>
#include <map>
#include <vector> #include <vector>
#include <dune/common/classname.hh>
#include <dune/common/exceptions.hh>
#include <dune/geometry/referenceelements.hh> #include <dune/geometry/referenceelements.hh>
#include <dune/pdelab/common/quadraturerules.hh> #include <dune/pdelab/common/quadraturerules.hh>
...@@ -39,6 +33,9 @@ namespace Operator { ...@@ -39,6 +33,9 @@ namespace Operator {
* \Gamma_D * \Gamma_D
* @f} * @f}
* *
* @todo Use diffusion coefficient grid function.
* @todo Implement outflow boundary condition.
*
* @tparam Boundary Type of the class providing boundary * @tparam Boundary Type of the class providing boundary
* conditions * conditions
* @tparam GridFunctionWaterFlux Type of a grid function which provides * @tparam GridFunctionWaterFlux Type of a grid function which provides
...@@ -182,20 +179,19 @@ public: ...@@ -182,20 +179,19 @@ public:
auto geo_i = entity_i.geometry(); auto geo_i = entity_i.geometry();
auto geo_o = entity_o.geometry(); auto geo_o = entity_o.geometry();
// Get geometry of intersection in local coordinates of neighbor entities // Get geometry of intersection in local coord. of neighbor entities
auto geo_in_i = ig.geometryInInside(); auto geo_in_i = ig.geometryInInside();
// Face volume for integration // Face volume for integration
auto face_volume = geo.volume(); auto face_volume = geo.volume();
// Evaluate diffusion coeff. from both sides and take harmonic average // Evaluate diffusion coeff. from both sides and take harmonic average
auto diff_coeff_i = _diff_coeff; // FIXME auto diff_coeff_i = _diff_coeff; // FIXME: evaluate from grid function
auto diff_coeff_o = _diff_coeff; // FIXME auto diff_coeff_o = _diff_coeff; // FIXME: evaluate from grid function
auto normal_face = ig.centerUnitOuterNormal(); auto normal_face = ig.centerUnitOuterNormal();
auto diff_coeff_avg = 2.0/( 1.0/(diff_coeff_i+ 1E-30) auto diff_coeff_avg = 2.0/( 1.0/(diff_coeff_i+ 1E-30)
+ 1.0/(diff_coeff_o+ 1E-30)); + 1.0/(diff_coeff_o+ 1E-30));
// Get center position of the face w.r.t inside entity // Get center position of the face w.r.t inside entity
auto p_center_face_i = geo_in_i.center(); auto p_center_face_i = geo_in_i.center();
...@@ -204,18 +200,18 @@ public: ...@@ -204,18 +200,18 @@ public:
_gf_flux->evaluate(entity_i, p_center_face_i, water_flux_i); _gf_flux->evaluate(entity_i, p_center_face_i, water_flux_i);
auto water_flux_n = water_flux_i*normal_face; auto water_flux_n = water_flux_i*normal_face;
// Inside solute value // Inside/outside solute value
RF u = x_i(lfsu_i,0); RF u_i = x_i(lfsu_i,0);
#if 1 RF u_o = x_o(lfsu_o,0);
// Upwinding // Upwinding
u = (water_flux_n>=0) ? x_i(lfsu_i,0) : x_o(lfsu_o,0); const auto& u_upwind = (water_flux_n>=0) ? u_i : u_o;
#endif
// Entity centers in references elements // Entity centers in references elements
auto ref_el_i = referenceElement(geo_i); const auto& ref_el_i = referenceElement(geo_i);
auto ref_el_o = referenceElement(geo_o); const auto& ref_el_o = referenceElement(geo_o);
auto p_center_i = ref_el_i.position(0,0); const auto& p_center_i = ref_el_i.position(0,0);
auto p_center_o = ref_el_o.position(0,0); const auto& p_center_o = ref_el_o.position(0,0);
// Entity centers in global coordinates // Entity centers in global coordinates
auto p_center_i_global = geo_i.global(p_center_i); auto p_center_i_global = geo_i.global(p_center_i);
...@@ -226,14 +222,14 @@ public: ...@@ -226,14 +222,14 @@ public:
auto distance = p_center_i_global.two_norm(); auto distance = p_center_i_global.two_norm();
// Finite difference of u between the two entities // Finite difference of u between the two entities
RF dudn = (x_o(lfsu_o,0)-x_i(lfsu_i,0))/distance; RF dudn = (u_o-u_i)/distance;
// Evaluate saturation // Evaluate saturation
typename GridFunctionSaturation::Traits::RangeType sat_i; typename GridFunctionSaturation::Traits::RangeType sat_i;
_gf_sat->evaluate(entity_i,p_center_i,sat_i); _gf_sat->evaluate(entity_i,p_center_i,sat_i);
// Solute flux in normal direction w.r.t the intersection // Solute flux in normal direction w.r.t the intersection
auto normal_flux = u*water_flux_n - diff_coeff_avg*dudn*sat_i; auto normal_flux = u_upwind*water_flux_n - diff_coeff_avg*dudn*sat_i;
// Symmetric contribution to residual on inside element // Symmetric contribution to residual on inside element
r_i.accumulate(lfsv_i, 0, normal_flux*face_volume ); r_i.accumulate(lfsv_i, 0, normal_flux*face_volume );
...@@ -267,8 +263,8 @@ public: ...@@ -267,8 +263,8 @@ public:
// Get geometries // Get geometries
auto geo = entitiy_face.geometry(); auto geo = entitiy_face.geometry();
auto ref_el = referenceElement(geo); const auto& ref_el = referenceElement(geo);
auto p_center_face = ref_el.position(0,0); const auto& p_center_face = ref_el.position(0,0);
auto bc = _boundary.bc(entitiy_face,p_center_face,_time); auto bc = _boundary.bc(entitiy_face,p_center_face,_time);
...@@ -309,7 +305,7 @@ public: ...@@ -309,7 +305,7 @@ public:
_gf_sat->evaluate(entity_i, p_center_face_i, sat_i); _gf_sat->evaluate(entity_i, p_center_face_i, sat_i);
// Evaluate diffusion coefficient from inside // Evaluate diffusion coefficient from inside
auto diff_coeff_i = _diff_coeff; auto diff_coeff_i = _diff_coeff; // FIXME: evaluate from grid function
// Inside unknown value // Inside unknown value
RF u = x_i(lfsu_i,0); RF u = x_i(lfsu_i,0);
...@@ -465,15 +461,15 @@ public: ...@@ -465,15 +461,15 @@ public:
auto geo = eg.geometry(); auto geo = eg.geometry();
// inside cell center // inside cell center
auto ref_el = referenceElement(geo); const auto& ref_el = referenceElement(geo);
auto p_center = ref_el.position(0,0); const auto& p_center = ref_el.position(0,0);
typename GridFunctionSaturation::Traits::RangeType saturation; typename GridFunctionSaturation::Traits::RangeType saturation;
_gfSaturation->evaluate(eg.entity(),p_center,saturation); _gfSaturation->evaluate(eg.entity(),p_center,saturation);
// update residual // update residual
r.accumulate(lfsu ,0 , saturation*u*geo.volume()); r.accumulate(lfsv ,0 , saturation*u*geo.volume());
} }
/*---------------------------------------------------------------------*//** /*---------------------------------------------------------------------*//**
......
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