From 7ba0f8fb104969146de12ed82d0b4caacca0d394 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 3 Jul 2017 18:05:04 +0200
Subject: [PATCH] [dt] Set current x to old x when a step is rejected
It affected the initial guess for linear solver in the Picard linearization.
---
NumLib/ODESolver/TimeDiscretization.h | 28 +++++++++++++++++++++++
ProcessLib/UncoupledProcessesTimeLoop.cpp | 7 +++---
2 files changed, 32 insertions(+), 3 deletions(-)
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index db3fb84d8b7..d14e2fd4180 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -142,6 +142,14 @@ public:
virtual void pushState(const double t, GlobalVector const& x,
InternalMatrixStorage const& strg) = 0;
+ /*! Indicate that the current timestep is rejected and that the computation
+ * of the present time step will be restarted with a new time step size.
+ *
+ * \param x The solution at the current time step, which is going to be
+ * reset to its previous value.
+ */
+ virtual void popState(GlobalVector& x) = 0;
+
/*! Indicate that the computation of a new timestep is being started now.
*
* \warning Currently changing timestep sizes are not supported. Thus,
@@ -262,6 +270,11 @@ public:
MathLib::LinAlg::copy(x, _x_old);
}
+ void popState(GlobalVector& x) override
+ {
+ MathLib::LinAlg::copy(_x_old, x);
+ }
+
void nextTimestep(const double t, const double delta_t) override
{
_t = t;
@@ -314,6 +327,11 @@ public:
MathLib::LinAlg::copy(x, _x_old);
}
+ void popState(GlobalVector& x) override
+ {
+ MathLib::LinAlg::copy(_x_old, x);
+ }
+
void nextTimestep(const double t, const double delta_t) override
{
_t_old = _t;
@@ -392,6 +410,11 @@ public:
strg.pushMatrices();
}
+ void popState(GlobalVector& x) override
+ {
+ MathLib::LinAlg::copy(_x_old, x);
+ }
+
void nextTimestep(const double t, const double delta_t) override
{
_t = t;
@@ -463,6 +486,11 @@ public:
void pushState(const double, GlobalVector const& x,
InternalMatrixStorage const&) override;
+ void popState(GlobalVector& x) override
+ {
+ MathLib::LinAlg::copy(*_xs_old[_xs_old.size()-1], x);
+ }
+
void nextTimestep(const double t, const double delta_t) override
{
_t = t;
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index 07ebb624fef..564f87ec24e 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -620,11 +620,11 @@ double UncoupledProcessesTimeLoop::computeTimeStepping(
continue;
}
+ auto& time_disc = ppd.time_disc;
+ auto& mat_strg = *ppd.mat_strg;
+ auto& x = *_process_solutions[i];
if (all_process_steps_accepted)
{
- auto& time_disc = ppd.time_disc;
- auto& mat_strg = *ppd.mat_strg;
- auto const& x = *_process_solutions[i];
time_disc->pushState(t, x, mat_strg);
}
else
@@ -635,6 +635,7 @@ double UncoupledProcessesTimeLoop::computeTimeStepping(
"Time step %d is rejected. "
"The computation is back to the previous time.",
accepted_steps + rejected_steps);
+ time_disc->popState(x);
}
}
}
--
GitLab