Commit 1517c757 by Lukas Riedel 📝 Committed by Santiago Ospina De Los Ríos

### Implement outflow BC in Richards model

* Outflow BC is a variant of the Dirichlet BC where only negative fluxes
(out of the domain) are applied. Positive fluxes are ignored and
effectively yield a no flow BC.
* Enable "Outflow" as Richards model BC.
parent 24d28d96
 ... @@ -37,6 +37,7 @@ ... @@ -37,6 +37,7 @@ * Control negative transport solution by a check policy !181 * Control negative transport solution by a check policy !181 * DG solver for solute transport model !112 * DG solver for solute transport model !112 * Cookbook tutorial on using the random field generator !184 * Cookbook tutorial on using the random field generator !184 * Outflow boundary condition for Richards model !191 ### Changed ### Changed * Data structures for storing and accessing parameter information !55 * Data structures for storing and accessing parameter information !55 ... ...
 ... @@ -98,85 +98,99 @@ A boundary condition has the following template: ... @@ -98,85 +98,99 @@ A boundary condition has the following template: value: value: # # .. option:: Dirichlet Dirichlet ^^^^^^^^^ A Dirichlet boundary conditions determines the exact value of the solution A Dirichlet boundary conditions determines the exact value of the solution at the boundary. at the boundary. **Richards Model** Richards Model Quantity: Matric head :math:h_m \, [\text{m}]. Quantity: Matric head :math:h_m \, [\text{m}]. .. tip:: .. tip:: Typical use cases for Dirichlet conditions in a Richards Model Typical use cases for Dirichlet conditions in a Richards Model include: include: * Height above a fixed water table at the lower boundary. * Height above a fixed water table at the lower boundary. * Evaporative potential at upper boundary, if :math:h_m < z_\text{wt}, * Evaporative potential at upper boundary, if :math:h_m < z_\text{wt}, where :math:z_\text{wt} is the height above the water table. where :math:z_\text{wt} is the height above the water table. **Transport Model** Transport Model Quantity: Total concentration :math:C_t\,[\text{kg}\,\text{m}^{-d}] Quantity: Total concentration :math:C_t\,[\text{kg}\,\text{m}^{-d}] or concentration in the water phase or concentration in the water phase :math:C_w\,[\text{kg}\,\text{m}^{-d}], depending on the value of :math:C_w\,[\text{kg}\,\text{m}^{-d}], depending on the value of concentration_type. The integer :math:d denotes the spatial concentration_type. The integer :math:d denotes the spatial dimension. dimension. Variables: Variables * type: Dirichlet * type: Dirichlet * value: Quantity as specified above. * value: Quantity as specified above. * concentration_type **(Transport only)**: Switch between total * concentration_type **(Transport only)**: Switch between total (:math:C_t) or water_phase (:math:C_w) concentration as (:math:C_t) or water_phase (:math:C_w) concentration as defined above. defined above. .. option:: Neumann Neumann ^^^^^^^ A Neumann boundary condition determines the value of target quantity flux A Neumann boundary condition determines the value of target quantity flux at the boundary. at the boundary. **Richards Model** Richards Model Quantity: Water flux :math:\mathbf{j}_w \, [\text{ms}^{-1}]. Quantity: Water flux :math:\mathbf{j}_w \, [\text{ms}^{-1}]. Precipitation and evapotranspiration fluxes are typically measured Precipitation and evapotranspiration fluxes are typically measured with respect to a flat surface. For modeling these fluxes, the with respect to a flat surface. For modeling these fluxes, the horizontal_projection key is used, which scales the flux applied horizontal_projection key is used, which scales the flux applied to the boundary with its inclination with respect to the direction of to the boundary with its inclination with respect to the direction of gravitational force, gravitational force, .. math:: .. math:: j_{N, \text{scaled}} j_{N, \text{scaled}} = \left| \mathbf{\hat{n}} \cdot \mathbf{\hat{g}} \right| j_N . = \left| \mathbf{\hat{n}} \cdot \mathbf{\hat{g}} \right| j_N . .. note:: Hydraulic conductivity decreases with water content. Evaporation .. note:: Hydraulic conductivity decreases with water content. Evaporation fluxes must not infinitely exceed the stable limit, or else the fluxes must not infinitely exceed the stable limit, or else the solver will inevitably fail. solver will inevitably fail. **Transport Model** Transport Model Quantity: Solute flux Quantity: Solute flux :math:\mathbf{j}_s \, [\text{kg}\,\text{m}^{-d+1}\,\text{s}^{-1}], :math:\mathbf{j}_s \, [\text{kg}\,\text{m}^{-d+1}\,\text{s}^{-1}], where :math:d indicates the spatial dimensions. where :math:d indicates the spatial dimensions. Variables * type: Neumann * value: Flux perpendicular to the surface :math:j_N, in the quantities defined above. Fluxes into the domain are negative, fluxes out of the domain positive. * horizontal_projection **(Richards only)**: Boolean if the actual boundary flux should be scaled with the surface inclination against the direction of gravity. Outflow ^^^^^^^ An outflow boundary condition only allows flow out of the domain. Variables: Richards Model * type: Neumann Quantity: Matric head :math:h_m \, [\text{m}] . * value: Flux perpendicular to the surface :math:j_N, in the quantities defined above. Fluxes into the domain are negative, fluxes out of the domain positive. * horizontal_projection **(Richards only)**: Boolean if the actual boundary flux should be scaled with the surface inclination against the direction of gravity. .. option:: Outflow (Transport only) The outflow boundary condition is a variant of the Dirichlet boundary condition. It applies the specified matric head only if it leads to a net flux out of the domain. Otherwise, the flux is set to zero (no flow boundary condition). An outflow boundary condition allows solute be transported throw the Transport Model Quantity: Solute flux :math:\mathbf{j}_s \, [\text{kg}\,\text{m}^{-d+1}\,\text{s}^{-1}], where :math:d indicates the spatial dimensions. The outflow boundary condition allows solute be transported throw the boundary by the water flux. As solute concentration and water flux are boundary by the water flux. As solute concentration and water flux are only modelled inside the domain, this effectively only allows the solute only modelled inside the domain, this effectively only allows the solute to leave the domain and be discarded. The boundary condition value is an to leave the domain and be discarded. The boundary condition value is an additional solute flux for which all considerations of a Neumann boundary additional solute flux for which all considerations of a Neumann condition apply (see above). boundary condition apply (see above). Quantity: Solute flux :math:\mathbf{j}_s \, [\text{kg}\,\text{m}^{-d+1}\,\text{s}^{-1}], where :math:d indicates the spatial dimensions. The total solute flux through the boundary is then given by The total solute flux through the boundary is then given by ... @@ -185,15 +199,16 @@ A boundary condition has the following template: ... @@ -185,15 +199,16 @@ A boundary condition has the following template: j_s = \left| \mathbf{j}_s \cdot \mathbf{\hat{n}} \right| + j_N , j_s = \left| \mathbf{j}_s \cdot \mathbf{\hat{n}} \right| + j_N , where :math:\mathbf{j}_s = C_w \mathbf{j}_w is the solute flux at the where :math:\mathbf{j}_s = C_w \mathbf{j}_w is the solute flux at the boundary, :math:\mathbf{\hat{n}} is the boundary unit normal vector and boundary, :math:\mathbf{\hat{n}} is the boundary unit normal vector :math:j_N is the given boundary condition value, as defined below. and :math:j_N is the flux perpendicular to the surface as specified in the boundary condition value. This flux is applied in addition to Variables: regular seepage through the boundary. It is recommended to set this * type: Outflow value to 0 (zero). * value: Flux perpendicular to the surface :math:j_N, as in :option:Neumann. This flux is applied in addition to regular Variables seepage through the boundary. It is recommended to set this value to * type: Outflow 0 (zero). * value: Matric head :math:h_m or solute flux perpendicular to the surface :math:j_N in quantities defined above, depending on the model. Example Files Example Files ------------- ------------- ... ...
 ... @@ -62,24 +62,11 @@ private: ... @@ -62,24 +62,11 @@ private: }; }; /// Return the type strings of permitted boundary conditions /// Return the type strings of permitted boundary conditions /** Currently, both models support all boundary condition types. */ static constexpr auto get_permitted_type_strs () static constexpr auto get_permitted_type_strs () { { if constexpr (mode == BCMMode::richards) { return known_type_strs; return std::array{ DirichletBoundaryCondition::type_str, NeumannBoundaryCondition::type_str }; } else if constexpr (mode == BCMMode::transport) { return std::array{ DirichletBoundaryCondition::type_str, NeumannBoundaryCondition::type_str, OutflowBoundaryCondition::type_str }; } else { return known_type_strs; } } } /// Return true of the given BC type is permitted for this model /// Return true of the given BC type is permitted for this model ... ...
 ... @@ -602,6 +602,7 @@ public: ... @@ -602,6 +602,7 @@ public: // query the boundary condition // query the boundary condition const auto bc = boundary->bc(ig.intersection()); const auto bc = boundary->bc(ig.intersection()); const auto bc_type = bc->type(); // loop over quadrature points and integrate normal flux // loop over quadrature points and integrate normal flux for (const auto& it : quadratureRule(gtface,intorder)) for (const auto& it : quadratureRule(gtface,intorder)) ... @@ -654,7 +655,7 @@ public: ... @@ -654,7 +655,7 @@ public: u_s += x_s(lfsu_s,i) * phiu_s[i]; u_s += x_s(lfsu_s,i) * phiu_s[i]; // update residual according to byType flag // update residual according to byType flag if (bc->type() == BCType::Neumann) if (bc_type == BCType::Neumann) { { // flux is given parallel to gravity vector // flux is given parallel to gravity vector RF normal_flux = bc->evaluate(time); RF normal_flux = bc->evaluate(time); ... @@ -673,7 +674,7 @@ public: ... @@ -673,7 +674,7 @@ public: continue; continue; } } else if (bc->type() == BCType::Dirichlet) else if (bc_type == BCType::Dirichlet or bc_type == BCType::Outflow) { { // compute gradient of u // compute gradient of u Vector gradu_s(0.); Vector gradu_s(0.); ... @@ -718,6 +719,11 @@ public: ... @@ -718,6 +719,11 @@ public: relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n); relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n); } } // Avoid accumulating inflow for outflow BC if (bc_type == BCType::Outflow and FloatCmp::lt(double(numFlux), 0.0)) { continue; } if constexpr (lopcase != LOPCase::RTVolume) if constexpr (lopcase != LOPCase::RTVolume) { { ... ...
 #ifndef DUNE_DORIE_RICHARDS_OPERATOR_FV_HH #ifndef DUNE_DORIE_RICHARDS_OPERATOR_FV_HH #define DUNE_DORIE_RICHARDS_OPERATOR_FV_HH #define DUNE_DORIE_RICHARDS_OPERATOR_FV_HH #include #include #include #include #include ... @@ -261,7 +263,8 @@ public: ... @@ -261,7 +263,8 @@ public: const auto saturation_f_i = _param->saturation_f(); const auto saturation_f_i = _param->saturation_f(); const auto conductivity_f_i = _param->conductivity_f(); const auto conductivity_f_i = _param->conductivity_f(); if (bc->type() == BCType::Dirichlet) const auto bc_type = bc->type(); if (bc_type == BCType::Dirichlet or bc_type == BCType::Outflow) { { const auto geo_i = entity_i.geometry(); const auto geo_i = entity_i.geometry(); const auto center_position_i_g = geo_i.center(); const auto center_position_i_g = geo_i.center(); ... @@ -291,10 +294,16 @@ public: ... @@ -291,10 +294,16 @@ public: // Water flux in normal direction w.r.t the intersection // Water flux in normal direction w.r.t the intersection const auto water_flux_n = - cond_i*dudn; const auto water_flux_n = - cond_i*dudn; // Avoid accumulation of inflow for outflow BC if (bc_type == BCType::Outflow and FloatCmp::lt(double(water_flux_n), 0.0)) { return; } // Contribution to residual from Dirichlet boundary // Contribution to residual from Dirichlet boundary r_i.accumulate(lfsv_i, 0, water_flux_n*volume_f); r_i.accumulate(lfsv_i, 0, water_flux_n*volume_f); } } else if (bc->type() == BCType::Neumann) else if (bc_type == BCType::Neumann) { { // Contribution to residual from Neumann boundary // Contribution to residual from Neumann boundary RangeFieldU water_flux_f = bc->evaluate(_time); RangeFieldU water_flux_f = bc->evaluate(_time); ... ...
 ... @@ -122,7 +122,7 @@ TEST_F(BCFactory, Outflow) ... @@ -122,7 +122,7 @@ TEST_F(BCFactory, Outflow) "time: 0.0" "time: 0.0" "}"); "}"); EXPECT_NO_THROW(_bcf_none->create("name", config, 1.0)); EXPECT_NO_THROW(_bcf_none->create("name", config, 1.0)); EXPECT_THROW(_bcf_richards->create("name", config, 1.0), Dune::IOError); EXPECT_NO_THROW(_bcf_richards->create("name", config, 1.0)); EXPECT_NO_THROW(_bcf_transport->create("name", config, 1.0)); EXPECT_NO_THROW(_bcf_transport->create("name", config, 1.0)); } } ... ...
 include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini _test_command = run _asset_path = "${CMAKE_CURRENT_LIST_DIR}" _evaluation = ode _order = 0, 1 | expand prec _upwind = none, semiUpwind, fullUpwind | expand _ic_mode = trans, stat | expand ic_mode __name = ode_outflow_{_order}_{_upwind}_{_ic_mode} [grid] gridType = rectangular initialLevel = 0 cells = 1 320, 1 160 | expand prec [richards] output.fileName = {__name} output.outputPath = {__name} output.vertexData = false initial.type = analytic, stationary | expand ic_mode initial.equation = -y boundary.file = "{_asset_path}/bcs/outflow.yml" parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E7, 0 | expand ic_mode time.maxTimestep = 1E7 time.startTimestep = 1E4 NewtonParameters.AbsoluteLimit = 1E-10 NewtonParameters.Reduction = 1E-10 numerics.penaltyFactor = 10 numerics.FEorder = {_order} numerics.upwinding = {_upwind} [_ode] headLower = 0. flux = -5.55e-6 head_abstol = 3E-4, 3.2E-5 | expand prec flux_abstol = 1E100, 3.1E-9 | expand prec flux_rt_abstol = 4E-14 # Error tolerance multipliers for stationary solution head_eps_stat = 1.0, 1.0 | expand prec flux_eps_stat = 1.0, 1.0 | expand prec flux_rt_eps_stat = 1E3