diff --git a/Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_initial_out_of_balance_forces.md b/Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_initial_out_of_balance_forces.md
new file mode 100644
index 0000000000000000000000000000000000000000..96a72b7a0b46e99d4abc8229254c406f4363ca1b
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_initial_out_of_balance_forces.md
@@ -0,0 +1 @@
+\copydoc NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_initial_out_of_balance_forces
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index c3216b1cc56db1714c0b3ab62148e3a73840d122..bf258e3d263e93eb4059f5a098fd5f3916121633 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -191,6 +191,17 @@ void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
// equation every time and could not forget it.
}
+void NonlinearSolver<NonlinearSolverTag::Newton>::calculateOutOfBalanceForces(
+ std::vector<GlobalVector*> const& x, int const process_id)
+{
+ if (!_compensate_initial_out_of_balance_forces)
+ return;
+
+ _equation_system->assemble(x, process_id);
+ _f_oob = &NumLib::GlobalVectorProvider::provider.getVector();
+ _equation_system->getResidual(*x[process_id], *_f_oob);
+}
+
NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const& postIterationCallback,
@@ -239,6 +250,10 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
sys.getJacobian(J);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());
+ // Apply out-of-balance forces if set.
+ if (_f_oob != nullptr)
+ LinAlg::axpy(res, -1, *_f_oob);
+
minus_delta_x.setZero();
timer_dirichlet.start();
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index bcffe0a4a476e2cbd6a2b108d461396ec2847c57..04c5815372afd60dc334dfca4541cd5f74718c9a 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -49,6 +49,9 @@ public:
virtual void assemble(std::vector<GlobalVector*> const& x,
int const process_id) const = 0;
+ virtual void calculateOutOfBalanceForces(
+ std::vector<GlobalVector*> const& x, int const process_id) = 0;
+
/*! Assemble and solve the equation system.
*
* \param x in: the initial guess, out: the solution.
@@ -113,12 +116,20 @@ public:
void assemble(std::vector<GlobalVector*> const& x,
int const process_id) const override;
+ void calculateOutOfBalanceForces(std::vector<GlobalVector*> const& x,
+ int const process_id) override;
+
NonlinearSolverStatus solve(
std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const&
postIterationCallback,
int const process_id) override;
+ void compensateInitialOutOfBalanceForces(bool const value)
+ {
+ _compensate_initial_out_of_balance_forces = value;
+ }
+
private:
GlobalLinearSolver& _linear_solver;
System* _equation_system = nullptr;
@@ -133,11 +144,18 @@ private:
//! conservative approach.
double const _damping;
+ GlobalVector* _f_oob = nullptr; //!< Out-of-balance forces vector.
std::size_t _res_id = 0u; //!< ID of the residual vector.
std::size_t _J_id = 0u; //!< ID of the Jacobian matrix.
std::size_t _minus_delta_x_id = 0u; //!< ID of the \f$ -\Delta x\f$ vector.
std::size_t _x_new_id =
0u; //!< ID of the vector storing \f$ x - (-\Delta x) \f$.
+
+ /// Enables computation of the out-of-balance forces \f$ F_{\rm oob} \f$
+ /// before the first time step. The forces are zero if the external forces
+ /// are in equilibrium with the initial state/initial conditions. During the
+ /// simulation the new residuum reads \f$ \tilde r = r - F_{\rm oob} \f$.
+ bool _compensate_initial_out_of_balance_forces = false;
};
/*! Find a solution to a nonlinear equation using the Picard fixpoint iteration
@@ -175,12 +193,25 @@ public:
void assemble(std::vector<GlobalVector*> const& x,
int const process_id) const override;
+ void calculateOutOfBalanceForces(GlobalVector const& /*x*/) override
+ {
+ if (!_compensate_initial_out_of_balance_forces)
+ return;
+
+ OGS_FATAL("Calculation of out-of-balance forces is not implemented.");
+ }
+
NonlinearSolverStatus solve(
std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const&
postIterationCallback,
int const process_id) override;
+ void compensateInitialOutOfBalanceForces(bool const value)
+ {
+ _compensate_initial_out_of_balance_forces = value;
+ }
+
private:
GlobalLinearSolver& _linear_solver;
System* _equation_system = nullptr;
@@ -193,6 +224,10 @@ private:
std::size_t _rhs_id = 0u; //!< ID of the right-hand side vector.
std::size_t _x_new_id = 0u; //!< ID of the vector storing the solution of
//! the linearized equation.
+
+ /// \copydoc
+ /// NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_initial_out_of_balance_forces
+ bool _compensate_initial_out_of_balance_forces = false;
};
/*! Creates a new nonlinear solver from the given configuration.
diff --git a/ProcessLib/CreateProcessData.cpp b/ProcessLib/CreateProcessData.cpp
index 4f46cae448f2593ebcf75a016966015269dabdbe..0c7955c5198b11fc056610e50ec5829ea5e30a52 100644
--- a/ProcessLib/CreateProcessData.cpp
+++ b/ProcessLib/CreateProcessData.cpp
@@ -23,7 +23,8 @@ static std::unique_ptr<ProcessData> makeProcessData(
int const process_id,
Process& process,
std::unique_ptr<NumLib::TimeDiscretization>&& time_disc,
- std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit)
+ std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit,
+ bool const compensate_initial_out_of_balance_forces)
{
using Tag = NumLib::NonlinearSolverTag;
@@ -31,6 +32,8 @@ static std::unique_ptr<ProcessData> makeProcessData(
dynamic_cast<NumLib::NonlinearSolver<Tag::Picard>*>(
&nonlinear_solver))
{
+ nonlinear_solver_picard->compensateInitialOutOfBalanceForces(
+ compensate_initial_out_of_balance_forces);
return std::make_unique<ProcessData>(
std::move(timestepper), *nonlinear_solver_picard,
std::move(conv_crit), std::move(time_disc), process_id, process);
@@ -39,6 +42,8 @@ static std::unique_ptr<ProcessData> makeProcessData(
dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
&nonlinear_solver))
{
+ nonlinear_solver_newton->compensateInitialOutOfBalanceForces(
+ compensate_initial_out_of_balance_forces);
return std::make_unique<ProcessData>(
std::move(timestepper), *nonlinear_solver_newton,
std::move(conv_crit), std::move(time_disc), process_id, process);
@@ -87,6 +92,11 @@ std::vector<std::unique_ptr<ProcessData>> createPerProcessData(
//! \ogs_file_param{prj__time_loop__processes__process__convergence_criterion}
pcs_config.getConfigSubtree("convergence_criterion"));
+ auto const compensate_initial_out_of_balance_forces =
+ //! \ogs_file_param{prj__time_loop__processes__process__compensate_initial_out_of_balance_forces}
+ pcs_config.getConfigParameter<bool>(
+ "compensate_initial_out_of_balance_forces", false);
+
//! \ogs_file_param{prj__time_loop__processes__process__output}
auto output = pcs_config.getConfigSubtreeOptional("output");
if (output)
@@ -103,7 +113,8 @@ std::vector<std::unique_ptr<ProcessData>> createPerProcessData(
per_process_data.emplace_back(
makeProcessData(std::move(timestepper), nl_slv, process_id, pcs,
- std::move(time_disc), std::move(conv_crit)));
+ std::move(time_disc), std::move(conv_crit),
+ compensate_initial_out_of_balance_forces));
++process_id;
}
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index e26104074fa8ed8264243d2668e216cc291ab6f5..af939682793e5554f421ef6ac39f376220bf5b2f 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -167,6 +167,31 @@ std::vector<GlobalVector*> setInitialConditions(
return process_solutions;
}
+void calculateOutOfBalanceForces(
+ std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
+ std::vector<GlobalVector*> process_solutions)
+{
+ INFO("Calculate out-of-balance forces.");
+ for (auto& process_data : per_process_data)
+ {
+ auto& ode_sys = *process_data->tdisc_ode_sys;
+
+ auto& time_disc = *process_data->time_disc;
+ auto& nonlinear_solver = process_data->nonlinear_solver;
+ auto& conv_crit = *process_data->conv_crit;
+
+ auto const nl_tag = process_data->nonlinear_solver_tag;
+ setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
+ // dummy values to handle the time derivative terms more or less
+ // correctly, i.e. to ignore them.
+ double const t = 0;
+ double const dt = 1;
+ time_disc.nextTimestep(t, dt);
+ nonlinear_solver.calculateOutOfBalanceForces(process_solutions,
+ process_data->process_id);
+ }
+}
+
NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
std::vector<GlobalVector*>& x, std::size_t const timestep, double const t,
double const delta_t, ProcessData const& process_data,
@@ -445,6 +470,7 @@ bool TimeLoop::loop()
bool const is_staggered_coupling =
!isMonolithicProcess(*_per_process_data[0]);
+ bool out_of_balance_forces_computed = false;
double t = _start_time;
std::size_t accepted_steps = 0;
std::size_t rejected_steps = 0;
@@ -472,6 +498,12 @@ bool TimeLoop::loop()
t, process_data->process_id);
}
+ if (!out_of_balance_forces_computed)
+ {
+ calculateOutOfBalanceForces(_per_process_data, _process_solutions);
+ out_of_balance_forces_computed = true;
+ }
+
if (is_staggered_coupling)
{
nonlinear_solver_status =