Implementation of nodal source terms

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/t_geometry.md

@@ -1 +1 @@
-The name of the on which the boundary condition is defined.
+The name of the geometry on which the boundary condition is defined.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/i_source_terms.md

@@ -0,0 +1 @@
+Contains source terms for the process variable.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Nodal/c_Nodal.md

@@ -0,0 +1 @@
+Declares a source term to be associated to a node.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Nodal/t_value.md

@@ -0,0 +1 @@
+The value for the nodal source term.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/i_source_term.md

@@ -0,0 +1 @@
+Defines a source term for a process variable.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_component.md

@@ -0,0 +1,3 @@
+For the multi-component process variables specifies the component of the
+source term. The component ids start with 0. For single-component
+variables this is optional and defaults to 0.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_geometrical_set.md

@@ -0,0 +1,4 @@
+The name of the geometrical set (see \ref ogs_file_param__gml__name) in which
+the \ref
+ogs_file_param__prj__process_variables__process_variable__source_terms__source_term__geometry
+"geometry" is defined.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_geometry.md

@@ -0,0 +1 @@
+The name of the geometry on which the source term is defined.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_type.md

@@ -0,0 +1 @@
+Sets the type of the source term.

MeshGeoToolsLib/CreateSearchLength.cpp

@@ -0,0 +1,49 @@
+/**
+ * @file
+ * @copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "CreateSearchLength.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Error.h"
+
+#include "MeshGeoToolsLib/HeuristicSearchLength.h"
+#include "MeshGeoToolsLib/SearchLength.h"
+
+namespace MeshGeoToolsLib
+{
+std::unique_ptr<MeshGeoToolsLib::SearchLength> createSearchLengthAlgorithm(
+    BaseLib::ConfigTree const& external_config, MeshLib::Mesh const& mesh)
+{
+    boost::optional<BaseLib::ConfigTree> config =
+        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm}
+        external_config.getConfigSubtreeOptional("search_length_algorithm");
+
+    if (!config)
+        return std::unique_ptr<MeshGeoToolsLib::SearchLength>{
+            new MeshGeoToolsLib::SearchLength()};
+
+    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm__type}
+    std::string const type = config->getConfigParameter<std::string>("type");
+
+    if (type == "fixed")
+    {
+        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm__fixed__value}
+        double const length = config->getConfigParameter<double>("value");
+        return std::unique_ptr<MeshGeoToolsLib::SearchLength>{
+            new MeshGeoToolsLib::SearchLength(length)};
+    }
+    if (type == "heuristic")
+    {
+        return std::unique_ptr<MeshGeoToolsLib::HeuristicSearchLength>{
+            new MeshGeoToolsLib::HeuristicSearchLength(mesh)};
+    }
+    OGS_FATAL("Unknown search length algorithm type '%s'.", type.c_str());
+}
+
+}  // end namespace MeshGeoToolsLib

MeshGeoToolsLib/CreateSearchLength.h

@@ -13,40 +13,20 @@
 
 #include <memory>
 
-#include "BaseLib/ConfigTree.h"
+namespace BaseLib
+{
+class ConfigTree;
+}
 
-#include "MeshGeoToolsLib/HeuristicSearchLength.h"
-#include "MeshGeoToolsLib/SearchLength.h"
+namespace MeshLib
+{
+class Mesh;
+}
 
 namespace MeshGeoToolsLib
 {
-std::unique_ptr<MeshGeoToolsLib::SearchLength> createSearchLengthAlgorithm(
-    BaseLib::ConfigTree const& external_config, MeshLib::Mesh const& mesh)
-{
-    boost::optional<BaseLib::ConfigTree> config =
-        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm}
-        external_config.getConfigSubtreeOptional("search_length_algorithm");
-
-    if (!config)
-        return std::unique_ptr<MeshGeoToolsLib::SearchLength>{
-            new MeshGeoToolsLib::SearchLength()};
-
-    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm__type}
-    std::string const type = config->getConfigParameter<std::string>("type");
-
-    if (type == "fixed")
-    {
-        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm__fixed__value}
-        double const length = config->getConfigParameter<double>("value");
-        return std::unique_ptr<MeshGeoToolsLib::SearchLength>{
-            new MeshGeoToolsLib::SearchLength(length)};
-    }
-    if (type == "heuristic")
-    {
-        return std::unique_ptr<MeshGeoToolsLib::HeuristicSearchLength>{
-            new MeshGeoToolsLib::HeuristicSearchLength(mesh)};
-    }
-    OGS_FATAL("Unknown search length algorithm type '%s'.", type.c_str());
-}
+class SearchLength;
 
+std::unique_ptr<MeshGeoToolsLib::SearchLength> createSearchLengthAlgorithm(
+    BaseLib::ConfigTree const& external_config, MeshLib::Mesh const& mesh);
 }  // end namespace MeshGeoToolsLib

ProcessLib/CMakeLists.txt

@@ -21,6 +21,7 @@ APPEND_SOURCE_FILES(SOURCES Parameter)
 APPEND_SOURCE_FILES(SOURCES PhaseField)
 APPEND_SOURCE_FILES(SOURCES RichardsFlow)
 APPEND_SOURCE_FILES(SOURCES SmallDeformation)
+APPEND_SOURCE_FILES(SOURCES SourceTerms)
 APPEND_SOURCE_FILES(SOURCES TES)
 APPEND_SOURCE_FILES(SOURCES ThermoMechanics)
 APPEND_SOURCE_FILES(SOURCES TwoPhaseFlowWithPP)

ProcessLib/GroundwaterFlow/Tests.cmake

@@ -597,3 +597,106 @@ AddTest(
     DIFF_DATA
     expected_dirichlet_nonuniform_pcs_0_ts_1_t_1.000000.vtu dirichlet_nonuniform_pcs_0_ts_1_t_1.000000.vtu pressure pressure
 )
+
+# tests for nodal source term implementation
+# For the special setup with a dirac source term at position (xi, eta) the
+# analytical solution in 2 dimensions is valid:
+# u(x,y) = ln(sqrt((x-xi)^2+(y-eta)^2))/(2 * Pi)
+AddTest(
+    NAME GroundWaterFlowProcess_NodalSourceTerm_circle_1e1
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e1_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 0.7 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e1_expected.vtu circle_1e1_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e1_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME GroundWaterFlowProcess_NodalSourceTerm_circle_1e2
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e2_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 1.1 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e2_expected.vtu circle_1e2_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e2_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME GroundWaterFlowProcess_NodalSourceTerm_circle_1e3
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e3_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 1.6 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e3_expected.vtu circle_1e3_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e3_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME GroundWaterFlowProcess_NodalSourceTerm_circle_1e4
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e4_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 1.8 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e4_expected.vtu circle_1e4_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e4_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME LARGE_GroundWaterFlowProcess_NodalSourceTerm_circle_1e5
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e5_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 2.15 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e5_expected.vtu circle_1e5_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e5_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME LARGE_GroundWaterFlowProcess_NodalSourceTerm_circle_1e6
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e6_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 2.52 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e6_expected.vtu circle_1e6_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e6_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME LARGE_GroundWaterFlowProcess_NodalSourceTerm_square_1e6
+    PATH Elliptic/square_1x1_GroundWaterFlow
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e6_with_nodal_sources.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 1.41 RELTOL 1e-16
+    DIFF_DATA
+    square_1x1_quad_1e6_nodal_sources_expected.vtu square_1e6_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS square_1e6_pcs_0_ts_1_t_1.000000.vtu
+)
+

ProcessLib/Process.cpp

@@ -59,6 +59,19 @@ void Process::initialize()
     DBUG("Initialize boundary conditions.");
     _boundary_conditions.addBCsForProcessVariables(
         _process_variables, *_local_to_global_index_map, _integration_order);
+
+    for (int variable_id = 0;
+         variable_id < static_cast<int>(_process_variables.size());
+         ++variable_id)
+    {
+        ProcessVariable& pv = _process_variables[variable_id];
+        auto sts =
+            pv.createSourceTerms(*_local_to_global_index_map, variable_id,
+                                 _integration_order);
+
+        std::move(sts.begin(), sts.end(),
+                  std::back_inserter(_source_terms));
+    }
 }
 
 void Process::setInitialConditions(double const t, GlobalVector& x)
@@ -128,6 +141,11 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
     assembleConcreteProcess(t, x, M, K, b);
 
     _boundary_conditions.applyNaturalBC(t, x, K, b);
+
+    for (auto const& st : _source_terms)
+    {
+        st->integrateNodalSourceTerm(t, b);
+    }
 }
 
 void Process::assembleWithJacobian(const double t, GlobalVector const& x,

ProcessLib/Process.h

@@ -14,6 +14,7 @@
 #include "NumLib/ODESolver/TimeDiscretization.h"
 #include "NumLib/NamedFunctionCaller.h"
 #include "ProcessLib/BoundaryCondition/BoundaryConditionCollection.h"
+#include "ProcessLib/SourceTerms/NodalSourceTerm.h"
 #include "ProcessLib/Parameter/Parameter.h"
 
 #include "ExtrapolatorData.h"
@@ -218,6 +219,7 @@ class Process
     std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
 
     BoundaryConditionCollection _boundary_conditions;
+    std::vector<std::unique_ptr<NodalSourceTerm>> _source_terms;
 
     ExtrapolatorData _extrapolator_data;
 };

ProcessLib/ProcessVariable.cpp

@@ -33,7 +33,8 @@ ProcessVariable::ProcessVariable(
           //! \ogs_file_param{prj__process_variables__process_variable__initial_condition}
           config.getConfigParameter<std::string>("initial_condition"),
           parameters, _n_components)),
-      _bc_builder(std::make_unique<BoundaryConditionBuilder>())
+      _bc_builder(std::make_unique<BoundaryConditionBuilder>()),
+      _source_term_builder(std::make_unique<SourceTermBuilder>())
 {
     DBUG("Constructing process variable %s", _name.c_str());
 
@@ -82,6 +83,49 @@ ProcessVariable::ProcessVariable(
     } else {
         INFO("No boundary conditions found.");
     }
+
+    // Source terms
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms}
+    if (auto sts_config = config.getConfigSubtreeOptional("source_terms"))
+    {
+        for (auto st_config :
+             //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term}
+             sts_config->getConfigSubtreeList("source_term"))
+        {
+            auto const geometrical_set_name =
+                    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__geometrical_set}
+                   st_config.getConfigParameter<std::string>("geometrical_set");
+            auto const geometry_name =
+                    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__geometry}
+                    st_config.getConfigParameter<std::string>("geometry");
+
+            GeoLib::GeoObject const* const geometry =
+                geometries.getGeoObject(geometrical_set_name, geometry_name);
+
+            if (! geometry)
+                OGS_FATAL(
+                    "No geometry with name `%s' has been found in the "
+                    "geometrical set `%s'.",
+                    geometry_name.c_str(), geometrical_set_name.c_str());
+
+            DBUG(
+                "Found geometry type \"%s\"",
+                GeoLib::convertGeoTypeToString(geometry->getGeoType()).c_str());
+
+            auto component_id =
+                //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__component}
+                st_config.getConfigParameterOptional<int>("component");
+
+            if (!component_id && _n_components == 1)
+                // default value for single component vars.
+                component_id = 0;
+
+            _source_term_configs.emplace_back(std::move(st_config), *geometry,
+                                              component_id);
+        }
+    } else {
+        INFO("No source terms found.");
+    }
 }
 
 ProcessVariable::ProcessVariable(ProcessVariable&& other)
@@ -91,7 +135,9 @@ ProcessVariable::ProcessVariable(ProcessVariable&& other)
       _shapefunction_order(other._shapefunction_order),
       _initial_condition(std::move(other._initial_condition)),
       _bc_configs(std::move(other._bc_configs)),
-      _bc_builder(std::move(other._bc_builder))
+      _bc_builder(std::move(other._bc_builder)),
+      _source_term_configs(std::move(other._source_term_configs)),
+      _source_term_builder(std::move(other._source_term_builder))
 {
 }
 
@@ -128,4 +174,20 @@ ProcessVariable::createBoundaryConditions(
     return bcs;
 }
 
+std::vector<std::unique_ptr<NodalSourceTerm>>
+ProcessVariable::createSourceTerms(
+    const NumLib::LocalToGlobalIndexMap& dof_table,
+    const int variable_id,
+    unsigned const integration_order)
+{
+    std::vector<std::unique_ptr<NodalSourceTerm>> source_terms;
+
+    for (auto& config : _source_term_configs)
+        source_terms.emplace_back(_source_term_builder->createSourceTerm(
+            config, dof_table, _mesh, variable_id, integration_order,
+            _shapefunction_order));
+
+    return source_terms;
+}
+
 }  // namespace ProcessLib