diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index 91ce5d379e0ccb152eb7f7e01a40025b69d7887c..8f2a5f25f43bb368be018544eac7ce0a87271691 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -73,7 +73,8 @@ void LocalAssemblerInterface::assembleWithJacobianForStaggeredScheme(
void LocalAssemblerInterface::computeSecondaryVariable(
std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const& dof_table, double const t,
- GlobalVector const& x, CoupledSolutionsForStaggeredScheme const* coupled_xs)
+ double const dt, GlobalVector const& x, GlobalVector const& x_dot,
+ CoupledSolutionsForStaggeredScheme const* coupled_xs)
{
auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
@@ -83,7 +84,8 @@ void LocalAssemblerInterface::computeSecondaryVariable(
}
auto const local_x = x.get(indices);
- computeSecondaryVariableConcrete(t, local_x);
+ auto const local_x_dot = x_dot.get(indices);
+ computeSecondaryVariableConcrete(t, dt, local_x, local_x_dot);
}
void LocalAssemblerInterface::setInitialConditions(
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index ef9a7c576f1cb947df1c6f5acef8d1379efa00bc..2b09c451092607e7990315e1b440c42d584a7d71 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -78,8 +78,9 @@ public:
virtual void computeSecondaryVariable(
std::size_t const mesh_item_id,
- NumLib::LocalToGlobalIndexMap const& dof_table, const double t,
- GlobalVector const& x,
+ NumLib::LocalToGlobalIndexMap const& dof_table, double const t,
+ double const dt, GlobalVector const& local_x,
+ GlobalVector const& local_x_dot,
CoupledSolutionsForStaggeredScheme const* coupled_xs);
virtual void preTimestep(std::size_t const mesh_item_id,
@@ -142,7 +143,10 @@ private:
}
virtual void computeSecondaryVariableConcrete(
- double const /*t*/, std::vector<double> const& /*local_x*/)
+ double const /*t*/,
+ double const /*dt*/,
+ std::vector<double> const& /*local_x*/,
+ std::vector<double> const& /*local_x_dot*/)
{
}
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 0c40e28cbfd343fbdb7fa04b2c942995ea0a291c..3fbc016c17bb49926c4cb21045231e5c000770f9 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -383,12 +383,16 @@ void Process::postNonLinearSolver(GlobalVector const& x, const double t,
postNonLinearSolverConcreteProcess(x, t, dt, process_id);
}
-void Process::computeSecondaryVariable(const double t, GlobalVector const& x,
+void Process::computeSecondaryVariable(double const t,
+ double const dt,
+ GlobalVector const& x,
+ GlobalVector const& x_dot,
int const process_id)
{
MathLib::LinAlg::setLocalAccessibleVector(x);
+ MathLib::LinAlg::setLocalAccessibleVector(x_dot);
- computeSecondaryVariableConcrete(t, x, process_id);
+ computeSecondaryVariableConcrete(t, dt, x, x_dot, process_id);
}
void Process::preIteration(const unsigned iter, const GlobalVector& x)
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index a7916c5daf9d4d4efbcbb99ebc5a195dcd1e0c99..5d5d92e32758fc7ccb849c12f2c1b1081c9b8d08 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -71,7 +71,10 @@ public:
void preIteration(const unsigned iter, GlobalVector const& x) final;
/// compute secondary variables for the coupled equations or for output.
- void computeSecondaryVariable(const double t, GlobalVector const& x,
+ void computeSecondaryVariable(double const t,
+ double const dt,
+ GlobalVector const& x,
+ GlobalVector const& x_dot,
int const process_id);
NumLib::IterationResult postIteration(GlobalVector const& x) final;
@@ -238,8 +241,10 @@ private:
{
}
- virtual void computeSecondaryVariableConcrete(const double /*t*/,
+ virtual void computeSecondaryVariableConcrete(double const /*t*/,
+ double const /*dt*/,
GlobalVector const& /*x*/,
+ GlobalVector const& /*x_dot*/,
int const /*process_id*/)
{
}
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 2c566463da162ada0b77bbc5630e7bdce0f7bb2b..8ab2d6f8f1e7d1a8400a616d1ffe22ee6303e58e 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -612,8 +612,25 @@ static constexpr std::string_view timestepper_cannot_reduce_dt =
void postTimestepForAllProcesses(
double const t, double const dt,
std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
- std::vector<GlobalVector*> const& process_solutions)
+ std::vector<GlobalVector*> const& process_solutions,
+ std::vector<GlobalVector*> const& process_solutions_prev)
{
+ std::vector<GlobalVector*> x_dots;
+ x_dots.reserve(per_process_data.size());
+ for (auto& process_data : per_process_data)
+ {
+ auto const process_id = process_data->process_id;
+ auto const& ode_sys = *process_data->tdisc_ode_sys;
+ auto const& time_discretization = *process_data->time_disc;
+
+ x_dots.emplace_back(&NumLib::GlobalVectorProvider::provider.getVector(
+ ode_sys.getMatrixSpecifications(process_id)));
+
+ time_discretization.getXdot(*process_solutions[process_id],
+ *process_solutions_prev[process_id],
+ *x_dots[process_id]);
+ }
+
// All _per_process_data share the first process.
bool const is_staggered_coupling =
!isMonolithicProcess(*per_process_data[0]);
@@ -630,8 +647,9 @@ void postTimestepForAllProcesses(
pcs.setCoupledSolutionsForStaggeredScheme(&coupled_solutions);
}
auto& x = *process_solutions[process_id];
+ auto& x_dot = *x_dots[process_id];
pcs.postTimestep(process_solutions, t, dt, process_id);
- pcs.computeSecondaryVariable(t, x, process_id);
+ pcs.computeSecondaryVariable(t, dt, x, x_dot, process_id);
}
}
@@ -668,7 +686,8 @@ NumLib::NonlinearSolverStatus TimeLoop::solveUncoupledEquationSystems(
}
}
- postTimestepForAllProcesses(t, dt, _per_process_data, _process_solutions);
+ postTimestepForAllProcesses(t, dt, _per_process_data, _process_solutions,
+ _process_solutions_prev);
return nonlinear_solver_status;
}
@@ -797,7 +816,8 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
}
- postTimestepForAllProcesses(t, dt, _per_process_data, _process_solutions);
+ postTimestepForAllProcesses(t, dt, _per_process_data, _process_solutions,
+ _process_solutions_prev);
return nonlinear_solver_status;
}
@@ -816,6 +836,7 @@ void TimeLoop::outputSolutions(bool const output_initial_condition,
{
auto const process_id = process_data->process_id;
auto& pcs = process_data->process;
+ auto const& ode_sys = *process_data->tdisc_ode_sys;
// If nonlinear solver diverged, the solution has already been
// saved.
if (!process_data->nonlinear_solver_status.error_norms_met)
@@ -827,12 +848,22 @@ void TimeLoop::outputSolutions(bool const output_initial_condition,
if (output_initial_condition)
{
- pcs.preTimestep(_process_solutions, _start_time,
- process_data->timestepper->getTimeStep().dt(),
- process_id);
+ // 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;
+ process_data->time_disc->nextTimestep(t, dt);
+
+ auto& x_dot = NumLib::GlobalVectorProvider::provider.getVector(
+ ode_sys.getMatrixSpecifications(process_id));
+ x_dot.setZero();
+
+ pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
// Update secondary variables, which might be uninitialized, before
// output.
- pcs.computeSecondaryVariable(_start_time, x, process_id);
+ pcs.computeSecondaryVariable(_start_time, dt, x, x_dot, process_id);
+
+ NumLib::GlobalVectorProvider::provider.releaseVector(x_dot);
}
if (is_staggered_coupling)
{