Merge branch 'fix-is-linear-output-off-by-one' into 'master'

Fixing off-by-one-timestep error in OGS's output when <linear> is used.

See merge request ogs/ogs!4749

ProcessLib/ComponentTransport/ComponentTransportProcess.cpp

@@ -219,32 +219,6 @@ void ComponentTransportProcess::assembleConcreteProcess(
     _asm_mat_cache.assemble(t, dt, x, x_prev, process_id, M, K, b, dof_tables,
                             _global_assembler, _local_assemblers,
                             pv.getActiveElementIDs());
-
-    BaseLib::RunTime time_residuum;
-    time_residuum.start();
-
-    if (_use_monolithic_scheme)
-    {
-        auto const residuum = computeResiduum(dt, *x[0], *x_prev[0], M, K, b);
-        for (std::size_t variable_id = 0; variable_id < _residua.size();
-             ++variable_id)
-        {
-            transformVariableFromGlobalVector(
-                residuum, variable_id, dof_tables[0], *_residua[variable_id],
-                std::negate<double>());
-        }
-    }
-    else
-    {
-        auto const residuum =
-            computeResiduum(dt, *x[process_id], *x_prev[process_id], M, K, b);
-        transformVariableFromGlobalVector(residuum, 0, dof_tables[process_id],
-                                          *_residua[process_id],
-                                          std::negate<double>());
-    }
-
-    INFO("[time] Computing residuum flow rates took {:g} s",
-         time_residuum.elapsed());
 }
 
 void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
@@ -478,5 +452,68 @@ void ComponentTransportProcess::postTimestepConcreteProcess(
                             pv.getActiveElementIDs());
 }
 
+void ComponentTransportProcess::preOutputConcreteProcess(
+    const double t,
+    double const dt,
+    std::vector<GlobalVector*> const& x,
+    std::vector<GlobalVector*> const& x_prev,
+    int const process_id)
+{
+    auto const matrix_specification = getMatrixSpecifications(process_id);
+
+    std::size_t matrix_id = 0u;
+    auto& M = NumLib::GlobalMatrixProvider::provider.getMatrix(
+        matrix_specification, matrix_id);
+    auto& K = NumLib::GlobalMatrixProvider::provider.getMatrix(
+        matrix_specification, matrix_id);
+    auto& b =
+        NumLib::GlobalVectorProvider::provider.getVector(matrix_specification);
+
+    M.setZero();
+    K.setZero();
+    b.setZero();
+
+    assembleConcreteProcess(t, dt, x, x_prev, process_id, M, K, b);
+
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
+    if (_use_monolithic_scheme)
+    {
+        dof_tables.push_back(std::ref(*_local_to_global_index_map));
+    }
+    else
+    {
+        std::generate_n(
+            std::back_inserter(dof_tables), _process_variables.size(),
+            [&]() { return std::ref(*_local_to_global_index_map); });
+    }
+
+    BaseLib::RunTime time_residuum;
+    time_residuum.start();
+
+    if (_use_monolithic_scheme)
+    {
+        auto const residuum = computeResiduum(dt, *x[0], *x_prev[0], M, K, b);
+        for (std::size_t variable_id = 0; variable_id < _residua.size();
+             ++variable_id)
+        {
+            transformVariableFromGlobalVector(
+                residuum, variable_id, dof_tables[0], *_residua[variable_id],
+                std::negate<double>());
+        }
+    }
+    else
+    {
+        auto const residuum =
+            computeResiduum(dt, *x[process_id], *x_prev[process_id], M, K, b);
+        transformVariableFromGlobalVector(residuum, 0, dof_tables[process_id],
+                                          *_residua[process_id],
+                                          std::negate<double>());
+    }
+
+    INFO("[time] Computing residuum flow rates took {:g} s",
+         time_residuum.elapsed());
+}
+
 }  // namespace ComponentTransport
 }  // namespace ProcessLib

ProcessLib/ComponentTransport/ComponentTransportProcess.h

@@ -169,6 +169,11 @@ class ComponentTransportProcess final : public Process
         GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
+    void preOutputConcreteProcess(const double t, double const dt,
+                                  std::vector<GlobalVector*> const& x,
+                                  std::vector<GlobalVector*> const& x_prev,
+                                  int const process_id) override;
+
     ComponentTransportProcessData _process_data;
 
     std::vector<std::unique_ptr<ComponentTransportLocalAssemblerInterface>>

ProcessLib/Output/Output.cpp

@@ -342,7 +342,7 @@ void Output::doOutput(Process const& process,
                       int const iteration,
                       std::vector<GlobalVector*> const& xs) const
 {
-    if (_output_data_specification.isOutputStep(timestep, t))
+    if (isOutputStep(timestep, t))
     {
         doOutputAlways(process, process_id, timestep, t, iteration, xs);
     }
@@ -361,7 +361,7 @@ void Output::doOutputLastTimestep(Process const& process,
                                   int const iteration,
                                   std::vector<GlobalVector*> const& xs) const
 {
-    if (!_output_data_specification.isOutputStep(timestep, t))
+    if (!isOutputStep(timestep, t))
     {
         doOutputAlways(process, process_id, timestep, t, iteration, xs);
     }
@@ -417,6 +417,11 @@ void Output::doOutputNonlinearIteration(
     INFO("[time] Output took {:g} s.", time_output.elapsed());
 }
 
+bool Output::isOutputStep(int const timestep, double const t) const
+{
+    return _output_data_specification.isOutputStep(timestep, t);
+}
+
 std::vector<std::string> Output::getFileNamesForOutput() const
 {
     auto construct_filename = ranges::views::transform(

ProcessLib/Output/Output.h

@@ -80,6 +80,9 @@ class Output
                                     const double t, const int iteration,
                                     std::vector<GlobalVector*> const& xs) const;
 
+    //! Tells if output will be written at the specified timestep/time.
+    bool isOutputStep(int const timestep, double const t) const;
+
     std::vector<double> const& getFixedOutputTimes() const
     {
         return _output_data_specification.fixed_output_times;

ProcessLib/Process.cpp

@@ -450,6 +450,17 @@ NumLib::IterationResult Process::postIteration(const GlobalVector& x)
     return postIterationConcreteProcess(x);
 }
 
+void Process::preOutput(const double t, double const dt,
+                        std::vector<GlobalVector*> const& x,
+                        std::vector<GlobalVector*> const& x_prev,
+                        int const process_id)
+{
+    for (auto const* solution : x)
+        MathLib::LinAlg::setLocalAccessibleVector(*solution);
+
+    preOutputConcreteProcess(t, dt, x, x_prev, process_id);
+}
+
 std::vector<GlobalIndexType>
 Process::getIndicesOfResiduumWithoutInitialCompensation() const
 {

ProcessLib/Process.h

@@ -136,6 +136,17 @@ class Process
                               GlobalMatrix& K, GlobalVector& b,
                               GlobalMatrix& Jac) final;
 
+    /* Computes data necessary for output.
+     *
+     * \pre Must be called before postTimestep() since the latter will overwrite
+     * previous states with current ones possibly affecting the data computed by
+     * this method.
+     */
+    void preOutput(const double t, double const dt,
+                   std::vector<GlobalVector*> const& x,
+                   std::vector<GlobalVector*> const& x_prev,
+                   int const process_id);
+
     std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
     getKnownSolutions(double const t, GlobalVector const& x,
                       int const process_id) const final
@@ -287,6 +298,13 @@ class Process
         return NumLib::IterationResult::SUCCESS;
     }
 
+    virtual void preOutputConcreteProcess(
+        const double /*t*/, double const /*dt*/,
+        std::vector<GlobalVector*> const& /*x*/,
+        std::vector<GlobalVector*> const& /*x_prev*/, int const /*process_id*/)
+    {
+    }
+
 protected:
     /** This function is for general cases, in which all equations of the
      coupled processes have the same number of unknowns. For the general cases

ProcessLib/TimeLoop.cpp

@@ -10,6 +10,8 @@
 
 #include "TimeLoop.h"
 
+#include <range/v3/algorithm/any_of.hpp>
+
 #include "BaseLib/Error.h"
 #include "BaseLib/RunTime.h"
 #include "CoupledSolutionsForStaggeredScheme.h"
@@ -38,6 +40,35 @@ void updateDeactivatedSubdomains(
     }
 }
 
+bool isOutputStep(std::vector<ProcessLib::Output> const& outputs,
+                  const int timestep, const double t)
+{
+    return ranges::any_of(outputs, [timestep, t](auto const& output)
+                          { return output.isOutputStep(timestep, t); });
+}
+
+void preOutputForAllProcesses(
+    int const timestep, double const t, double const dt,
+    std::vector<std::unique_ptr<ProcessLib::ProcessData>> const&
+        per_process_data,
+    std::vector<GlobalVector*> const& process_solutions,
+    std::vector<GlobalVector*> const& process_solutions_prev,
+    std::vector<ProcessLib::Output> const& outputs)
+{
+    if (!isOutputStep(outputs, timestep, t))
+    {
+        return;
+    }
+
+    for (auto& process_data : per_process_data)
+    {
+        auto const process_id = process_data->process_id;
+        auto& pcs = process_data->process;
+
+        pcs.preOutput(t, dt, process_solutions, process_solutions_prev,
+                      process_id);
+    }
+}
 }  // namespace
 
 namespace ProcessLib
@@ -217,6 +248,9 @@ NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
     auto const post_iteration_callback =
         [&](int iteration, std::vector<GlobalVector*> const& x)
     {
+        // Note: We don't call the postNonLinearSolver(), preOutput(),
+        // computeSecondaryVariable() and postTimestep() hooks here. This might
+        // lead to some inconsistencies in the data compared to regular output.
         for (auto const& output : outputs)
         {
             output.doOutputNonlinearIteration(process, process_id, timestep, t,
@@ -592,6 +626,13 @@ bool TimeLoop::preTsNonlinearSolvePostTs(double const t, double const dt,
     // iteration, an exception thrown in assembly, for example.
     if (nonlinear_solver_status.error_norms_met)
     {
+        // Later on, the timestep_algorithm might reject the timestep. We assume
+        // that this is a rare case, so still, we call preOutput() here. We
+        // don't expect a large overhead from it.
+        preOutputForAllProcesses(timesteps, t, dt, _per_process_data,
+                                 _process_solutions, _process_solutions_prev,
+                                 _outputs);
+
         postTimestepForAllProcesses(t, dt, _per_process_data,
                                     _process_solutions,
                                     _process_solutions_prev);
@@ -876,6 +917,10 @@ void TimeLoop::preOutputInitialConditions(const double t) const
             process_data->time_disc->nextTimestep(t, dt);
 
             pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
+
+            pcs.preOutput(_start_time, dt, _process_solutions,
+                          _process_solutions_prev, process_id);
+
             // Update secondary variables, which might be uninitialized, before
             // output.
             pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,
@@ -899,6 +944,10 @@ void TimeLoop::preOutputInitialConditions(const double t) const
             process_data->time_disc->nextTimestep(t, dt);
 
             pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
+
+            pcs.preOutput(_start_time, dt, _process_solutions,
+                          _process_solutions_prev, process_id);
+
             // Update secondary variables, which might be uninitialized, before
             // output.
             pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,

Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/DecayChain/GlobalImplicitApproach/1d_decay_chain_GIA.prj

@@ -240,6 +240,7 @@
                 <variable>pressure</variable>
                 <variable>LiquidMassFlowRate</variable>
                 <variable>[Cm-247]FlowRate</variable>
+                <variable>[Am-243]FlowRate</variable>
                 <variable>[Pu-239]FlowRate</variable>
                 <variable>[U-235]FlowRate</variable>
                 <variable>[Pa-231]FlowRate</variable>
@@ -497,7 +498,7 @@
         </vtkdiff>
         <vtkdiff>
             <regex>1d_decay_chain_GIA_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
-            <field>[Cm-247]FlowRate</field>
+            <field>[Am-243]FlowRate</field>
             <absolute_tolerance>1e-10</absolute_tolerance>
             <relative_tolerance>1e-16</relative_tolerance>
         </vtkdiff>