Merge branch 'NumLibChangeDirichletBcSettingForNewton' into 'master'

NumLib: Update setting of DirichletBCs in Newton See merge request ogs/ogs!4912
ufz · Jul 3, 2024 · 129f5a4 · 129f5a4
2 parents 6b3551e + d28e071
commit 129f5a4
Show file tree

Hide file tree

Showing 89 changed files with 874 additions and 885 deletions.
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -345,7 +345,6 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
 
         timer_dirichlet.start();
         sys.computeKnownSolutions(*x[process_id], process_id);
-        sys.applyKnownSolutions(*x[process_id]);
         time_dirichlet += timer_dirichlet.elapsed();
 
         sys.preIteration(iteration, *x[process_id]);
@@ -375,7 +374,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
         minus_delta_x.setZero();
 
         timer_dirichlet.start();
-        sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
+        sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
         time_dirichlet += timer_dirichlet.elapsed();
         INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
 

diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
@@ -75,9 +75,16 @@ class NonlinearSystem<NonlinearSolverTag::Newton> : public EquationSystem
     //! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
     //! \pre computeKnownSolutions() must have been called before.
     virtual void applyKnownSolutionsNewton(
-        GlobalMatrix& Jac, GlobalVector& res,
+        GlobalMatrix& Jac, GlobalVector& res, GlobalVector const& x,
         GlobalVector& minus_delta_x) const = 0;
 
+    //! Apply known solutions to the linearized equation system
+    //! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
+    //! \pre computeKnownSolutions() must have been called before.
+    virtual void applyKnownSolutionsPETScSNES(GlobalMatrix& Jac,
+                                              GlobalVector& res,
+                                              GlobalVector& x) const = 0;
+
     virtual void updateConstraints(GlobalVector& /*lower*/,
                                    GlobalVector& /*upper*/,
                                    int /*process_id*/) = 0;

diff --git a/NumLib/ODESolver/PETScNonlinearSolver.cpp b/NumLib/ODESolver/PETScNonlinearSolver.cpp
@@ -67,7 +67,7 @@ PetscErrorCode updateResidual(SNES /*snes*/, Vec x, Vec petsc_r,
     context->J->finalizeAssembly();
 
     context->system->getJacobian(*context->J);
-    context->system->applyKnownSolutionsNewton(
+    context->system->applyKnownSolutionsPETScSNES(
         *context->J, *context->r, *context->x[context->process_id]);
 
     VecCopy(context->r->getRawVector(), petsc_r);

diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -10,6 +10,9 @@
 
 #include "TimeDiscretizedODESystem.h"
 
+#include <range/v3/numeric/accumulate.hpp>
+#include <range/v3/view/transform.hpp>
+
 #include "MathLib/LinAlg/ApplyKnownSolution.h"
 #include "MathLib/LinAlg/UnifiedMatrixSetters.h"
 #include "NumLib/Exceptions.h"
@@ -133,13 +136,49 @@ void TimeDiscretizedODESystem<
 void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                               NonlinearSolverTag::Newton>::
     applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
+                              GlobalVector const& x,
                               GlobalVector& minus_delta_x) const
 {
     if (!_known_solutions)
     {
         return;
     }
 
+    using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+    std::size_t const size = ranges::accumulate(
+        *_known_solutions | ranges::views::transform([](auto const& bc)
+                                                     { return bc.ids.size(); }),
+        0);
+    std::vector<IndexType> ids;
+    ids.reserve(size);
+    std::vector<double> values;
+    values.reserve(size);
+
+    for (auto const& bc : *_known_solutions)
+    {
+        for (std::size_t i = 0; i < bc.ids.size(); ++i)
+        {
+            auto const id = bc.ids[i];
+            ids.push_back(id);
+            // minus_delta_x will be set to the difference between the current
+            // value and the Dirichlet BC value.
+            values.push_back(x[id] - bc.values[i]);
+        }
+    }
+
+    MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
+}
+
+void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
+                              NonlinearSolverTag::Newton>::
+    applyKnownSolutionsPETScSNES(GlobalMatrix& Jac, GlobalVector& res,
+                                 GlobalVector& x) const
+{
+    if (!_known_solutions)
+    {
+        return;
+    }
+
     using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
     std::vector<IndexType> ids;
     for (auto const& bc : *_known_solutions)
@@ -149,7 +188,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
 
     // For the Newton method the values must be zero
     std::vector<double> values(ids.size(), 0);
-    MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
+    MathLib::applyKnownSolution(Jac, res, x, ids, values);
 }
 
 TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,

diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -106,8 +106,12 @@ class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     void applyKnownSolutions(GlobalVector& x) const override;
 
     void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
+                                   GlobalVector const& x,
                                    GlobalVector& minus_delta_x) const override;
 
+    void applyKnownSolutionsPETScSNES(GlobalMatrix& Jac, GlobalVector& res,
+                                      GlobalVector& x) const override;
+
     void updateConstraints(GlobalVector& lower, GlobalVector& upper,
                            int const process_id) override
     {

diff --git a/ProcessLib/HydroMechanics/Tests.cmake b/ProcessLib/HydroMechanics/Tests.cmake
@@ -55,7 +55,7 @@ AddTest(
     TESTER vtkdiff
     REQUIREMENTS NOT OGS_USE_MPI
     DIFF_DATA
-    simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu pressure_interpolated pressure_interpolated 1e-10 0
+    simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu pressure_interpolated pressure_interpolated 1.1e-10 0
     simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu displacement displacement 1e-14 0
     simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu sigma sigma 1e-7 0
     simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu DarcyVelocity DarcyVelocity 1e-7 0
@@ -71,7 +71,7 @@ AddTest(
     TESTER vtkdiff
     REQUIREMENTS NOT OGS_USE_MPI
     DIFF_DATA
-    simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_python_ts_10_t_1000000.000000.vtu pressure_interpolated pressure_interpolated 1e-10 0
+    simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_python_ts_10_t_1000000.000000.vtu pressure_interpolated pressure_interpolated 1.1e-10 0
     simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_python_ts_10_t_1000000.000000.vtu displacement displacement 1e-14 0
     simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_python_ts_10_t_1000000.000000.vtu sigma sigma 1e-7 0
     simHM_ground_quadBCu_ts_10_t_1000000.000000.vtu simHM_ground_quadBCu_python_ts_10_t_1000000.000000.vtu DarcyVelocity DarcyVelocity 1e-7 0

diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/Tests.cmake b/ProcessLib/ThermalTwoPhaseFlowWithPP/Tests.cmake
@@ -1,11 +1,14 @@
+# Disabled until issue https://gitlab.opengeosys.org/ogs/ogs/-/issues/3486 is
+# resolved.
+set(ISSUE_3486_RESOLVED FALSE)
 AddTest(
     NAME 2D_Thermal_TwoPhase_heatpipe_small
     PATH Parabolic/ThermalTwoPhaseFlowPP/HeatPipe
     RUNTIME 10
     EXECUTABLE ogs
     EXECUTABLE_ARGS Twophase_HeatPipe_quad_curve_small.prj
     TESTER vtkdiff
-    REQUIREMENTS NOT OGS_USE_MPI
+    REQUIREMENTS NOT OGS_USE_MPI AND ISSUE_3486_RESOLVED
     DIFF_DATA
     ref_t_10000.000000.vtu thermaltwophaseflow_small_ts_100_t_10000.000000.vtu capillary_pressure capillary_pressure 1e-8 1e-10
     ref_t_10000.000000.vtu thermaltwophaseflow_small_ts_100_t_10000.000000.vtu gas_pressure gas_pressure 1e-8 1e-10
@@ -19,7 +22,7 @@ AddTest(
     EXECUTABLE ogs
     EXECUTABLE_ARGS Twophase_HeatPipe_quad_curve_large.prj
     TESTER vtkdiff
-    REQUIREMENTS NOT OGS_USE_MPI
+    REQUIREMENTS NOT OGS_USE_MPI AND ISSUE_3486_RESOLVED
     DIFF_DATA
     ref_t_1400000.000000.vtu thermaltwophaseflow_large_ts_1400_t_1400000.000000.vtu capillary_pressure capillary_pressure 9e-5 5e-10
     ref_t_1400000.000000.vtu thermaltwophaseflow_large_ts_1400_t_1400000.000000.vtu gas_pressure gas_pressure 5e-7 1e-10

diff --git a/ProcessLib/ThermoMechanics/Tests.cmake b/ProcessLib/ThermoMechanics/Tests.cmake
@@ -112,7 +112,7 @@ AddTest(
     REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
     RUNTIME 17
     DIFF_DATA
-    tm1_3Dgravity_ts_1_t_1.000000.vtu tm1_3Dgravity_ts_1_t_1.000000.vtu temperature temperature 1e-11 0.0
+    tm1_3Dgravity_ts_1_t_1.000000.vtu tm1_3Dgravity_ts_1_t_1.000000.vtu temperature temperature 4e-11 0.0
     tm1_3Dgravity_ts_1_t_1.000000.vtu tm1_3Dgravity_ts_1_t_1.000000.vtu displacement displacement 1e-7 0.0
 )
 

diff --git a/ProcessLib/ThermoRichardsFlow/Tests.cmake b/ProcessLib/ThermoRichardsFlow/Tests.cmake
@@ -149,7 +149,7 @@ AddTest(
     expected_TRMuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu temperature temperature 5e-5 1e-10
     expected_TRMuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu pressure pressure 5e-2 1e-6
     expected_TRMuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu saturation saturation 5e-5 1e-10
-    expected_TRMuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu MassFlowRate MassFlowRate 1e-9 1e-4
+    expected_TRMuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRuni_unsaturated_bishopstest_ts_10_t_1.000000.vtu MassFlowRate MassFlowRate 8e-6 1e-4
 )
 AddTest(
     NAME ThermoRichardsFlow_comp_TRMhyd_bishopstest-TRiso_bishopstest
@@ -164,7 +164,7 @@ AddTest(
     expected_TRMhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu temperature temperature 5e-5 1e-10
     expected_TRMhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu pressure pressure 5e-2 1e-6
     expected_TRMhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu saturation saturation 5e-5 1e-10
-    expected_TRMhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu MassFlowRate MassFlowRate 1e-9 1e-4
+    expected_TRMhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu TRhyd_unsaturated_bishopstest_ts_10_t_1.000000.vtu MassFlowRate MassFlowRate 5e-6 1e-4
 )
 AddTest(
     NAME ThermoRichardsFlow_TaskCDECOVALEX2023

diff --git a/Tests/Data/HydroMechanics/ExcavationNiches/excavation_niches.prj b/Tests/Data/HydroMechanics/ExcavationNiches/excavation_niches.prj
@@ -293,13 +293,13 @@
         <vtkdiff>
             <regex>excavation_niches_ts_.*.vtu</regex>
             <field>sigma</field>
-            <absolute_tolerance>4e-8</absolute_tolerance>
+            <absolute_tolerance>5e-8</absolute_tolerance>
             <relative_tolerance>1e-15</relative_tolerance>
         </vtkdiff>
         <vtkdiff>
             <regex>excavation_niches_ts_.*.vtu</regex>
             <field>pressure</field>
-            <absolute_tolerance>7e-8</absolute_tolerance>
+            <absolute_tolerance>8e-8</absolute_tolerance>
             <relative_tolerance>1e-15</relative_tolerance>
         </vtkdiff>
         <vtkdiff>

diff --git a/Tests/Data/LIE/Mechanics/coulomb_load_path_ts_0_t_0.000000.vtu b/Tests/Data/LIE/Mechanics/coulomb_load_path_ts_0_t_0.000000.vtu
@@ -30,6 +30,6 @@
     </Piece>
   </UnstructuredGrid>
   <AppendedData encoding="base64">
-   _AQAAAAAAAAAAgAAAAAAAABkAAAAAAAAAIQAAAAAAAAA=eF4z0zPVM9A1NjPSTTdKTrJIMzZN1UvJLCqpBABN8gcQAQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAADAAAAAAAAAA=eF5jYKAtAAAAYAABAQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAADAAAAAAAAAA=eF5jYKAtAAAAYAABAQAAAAAAAAAAgAAAAAAAAMAAAAAAAAAADAAAAAAAAAA=eF5jYBjaAAAAwAABAQAAAAAAAAAAgAAAAAAAADAAAAAAAAAAFgAAAAAAAAA=eF5jYMAPPmY87tj644kdjAYAYQALJw==AQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAADAAAAAAAAAA=eF5jYKAtAAAAYAABAQAAAAAAAAAAgAAAAAAAAMAAAAAAAAAADAAAAAAAAAA=eF5jYBjaAAAAwAABAQAAAAAAAAAAgAAAAAAAABgAAAAAAAAAEwAAAAAAAAA=eF77mPG4Y+uPJ3YMaAAAc2QFlA==AQAAAAAAAAAAgAAAAAAAADAAAAAAAAAADAAAAAAAAAA=eF5jYCANAAAAMAABAQAAAAAAAAAAgAAAAAAAADAAAAAAAAAADAAAAAAAAAA=eF5jYCANAAAAMAABAQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAADAAAAAAAAAA=eF5jYKAtAAAAYAABAQAAAAAAAAAAgAAAAAAAAJAAAAAAAAAAGwAAAAAAAAA=eF5jYMAHPtgTEHfALo9LHN08GB8hDgBcugZsAQAAAAAAAAAAgAAAAAAAAFAAAAAAAAAAHAAAAAAAAAA=eF5jZYAAFgZUwAilWdHkmdH4MHkmKA0ABegAIg==AQAAAAAAAAAAgAAAAAAAABgAAAAAAAAAEAAAAAAAAAA=eF5jYoAANijNBaUBAPgAEw==AQAAAAAAAAAAgAAAAAAAAAMAAAAAAAAACwAAAAAAAAA=eF5j5uQEAAAnABY=
+   _AQAAAAAAAAAAgAAAAAAAABkAAAAAAAAAIQAAAAAAAAA=eF4z0zPVM9A1NTDVTTe0tExONE9M1kvJLCqpBABM4AcPAQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAADAAAAAAAAAA=eF5jYKAtAAAAYAABAQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAADAAAAAAAAAA=eF5jYKAtAAAAYAABAQAAAAAAAAAAgAAAAAAAAMAAAAAAAAAADAAAAAAAAAA=eF5jYBjaAAAAwAABAQAAAAAAAAAAgAAAAAAAADAAAAAAAAAAFgAAAAAAAAA=eF5jYMAPPmY87tj644kdjAYAYQALJw==AQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAADAAAAAAAAAA=eF5jYKAtAAAAYAABAQAAAAAAAAAAgAAAAAAAAMAAAAAAAAAADAAAAAAAAAA=eF5jYBjaAAAAwAABAQAAAAAAAAAAgAAAAAAAABgAAAAAAAAAEwAAAAAAAAA=eF77mPG4Y+uPJ3YMaAAAc2QFlA==AQAAAAAAAAAAgAAAAAAAADAAAAAAAAAADAAAAAAAAAA=eF5jYCANAAAAMAABAQAAAAAAAAAAgAAAAAAAADAAAAAAAAAADAAAAAAAAAA=eF5jYCANAAAAMAABAQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAADAAAAAAAAAA=eF5jYKAtAAAAYAABAQAAAAAAAAAAgAAAAAAAAJAAAAAAAAAAGwAAAAAAAAA=eF5jYMAHPtgTEHfALo9LHN08GB8hDgBcugZsAQAAAAAAAAAAgAAAAAAAAFAAAAAAAAAAHAAAAAAAAAA=eF5jZYAAFgZUwAilWdHkmdH4MHkmKA0ABegAIg==AQAAAAAAAAAAgAAAAAAAABgAAAAAAAAAEAAAAAAAAAA=eF5jYoAANijNBaUBAPgAEw==AQAAAAAAAAAAgAAAAAAAAAMAAAAAAAAACwAAAAAAAAA=eF5j5uQEAAAnABY=
   </AppendedData>
 </VTKFile>