Commit b06aade2 authored by renchao.lu's avatar renchao.lu Committed by renchao.lu

wip.

parent efb2b965
/**
* \file
*
* \copyright
* Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#pragma once
#include "ControlNode.h"
#include <vector>
#include "BaseLib/RunTime.h"
#include "Output/Output.h"
namespace
{
void setEquationSystem(NumLib::NonlinearSolverBase& nonlinear_solver,
NumLib::EquationSystem& eq_sys,
NumLib::ConvergenceCriterion& conv_crit,
NumLib::NonlinearSolverTag nl_tag)
{
using Tag = NumLib::NonlinearSolverTag;
switch (nl_tag)
{
case Tag::Picard:
{
using EqSys = NumLib::NonlinearSystem<Tag::Picard>;
auto& eq_sys_ = static_cast<EqSys&>(eq_sys);
if (auto* nl_solver =
dynamic_cast<NumLib::NonlinearSolver<Tag::Picard>*>(
&nonlinear_solver);
nl_solver != nullptr)
{
nl_solver->setEquationSystem(eq_sys_, conv_crit);
}
else
{
OGS_FATAL(
"Could not cast nonlinear solver to Picard type solver.");
}
break;
}
case Tag::Newton:
{
using EqSys = NumLib::NonlinearSystem<Tag::Newton>;
auto& eq_sys_ = static_cast<EqSys&>(eq_sys);
if (auto* nl_solver =
dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
&nonlinear_solver);
nl_solver != nullptr)
{
nl_solver->setEquationSystem(eq_sys_, conv_crit);
}
#ifdef USE_PETSC
else if (auto* nl_solver =
dynamic_cast<NumLib::PETScNonlinearSolver*>(
&nonlinear_solver);
nl_solver != nullptr)
{
nl_solver->setEquationSystem(eq_sys_, conv_crit);
}
#endif // USE_PETSC
else
{
OGS_FATAL(
"Could not cast nonlinear solver to Newton type solver.");
}
break;
}
}
}
} // namespace
namespace ProcessLib
{
NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
std::vector<GlobalVector*>& x, std::vector<GlobalVector*> const& x_prev,
std::size_t const timestep, double const t, double const delta_t,
ProcessData const& process_data, Output& output_control)
{
auto& process = process_data.process;
int const process_id = process_data.process_id;
auto& time_disc = *process_data.time_disc;
auto& conv_crit = *process_data.conv_crit;
auto& ode_sys = *process_data.tdisc_ode_sys;
auto& nonlinear_solver = process_data.nonlinear_solver;
auto const nl_tag = process_data.nonlinear_solver_tag;
setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
// Note: Order matters!
// First advance to the next timestep, then set known solutions at that
// time, afterwards pass the right solution vector and time to the
// preTimestep() hook.
time_disc.nextTimestep(t, delta_t);
auto const post_iteration_callback =
[&](int iteration, std::vector<GlobalVector*> const& x) {
output_control.doOutputNonlinearIteration(
process, process_id, timestep, t, x, iteration);
};
auto const nonlinear_solver_status =
nonlinear_solver.solve(x, x_prev, post_iteration_callback, process_id);
if (nonlinear_solver_status.error_norms_met)
{
process.postNonLinearSolver(*x[process_id], t, delta_t, process_id);
}
return nonlinear_solver_status;
}
static NumLib::NonlinearSolverStatus solveMonolithicProcess(
const double t, const double dt, const std::size_t timestep_id,
ProcessData const& process_data, std::vector<GlobalVector*>& x,
std::vector<GlobalVector*> const& x_prev, Output& output)
{
BaseLib::RunTime time_timestep_process;
// time_timestep_process.start();
// auto const nonlinear_solver_status = solveOneTimeStepOneProcess(
// x, x_prev, timestep_id, t, dt, process_data, output);
// INFO("[time] Solving process #{:d} took {:g} s in time step #{:d} ",
// process_data.process_id, time_timestep_process.elapsed(),
// timestep_id);
// return nonlinear_solver_status;
}
NumLib::NonlinearSolverStatus ControlNode::solveEquationSystems(
const double t,
const double dt,
const std::size_t timestep_id,
std::unique_ptr<Output>& output)
{
NumLib::NonlinearSolverStatus nonlinear_solver_status;
bool const is_staggered_coupling = !isMonolithicProcess();
if (is_staggered_coupling)
{
nonlinear_solver_status = solveCoupledEquationSystemsByStaggeredScheme(
t, dt, timestep_id, output);
}
else
{
nonlinear_solver_status =
solveUncoupledEquationSystems(t, dt, timestep_id, output);
}
}
void ControlNode::calculateNonEquilibriumInitialResiduum()
{
INFO("Calculate non-equilibrium initial residuum.");
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.calculateNonEquilibriumInitialResiduum(
process_solutions, process_solutions_prev,
process_data->process_id);
}
}
void ControlNode::setCoupledSolutions()
{
bool const is_staggered_coupling = !isMonolithicProcess();
if (is_staggered_coupling)
{
for (auto& process_data : per_process_data)
{
auto const& x = *process_solutions[process_data->process_id];
// Create a vector to store the solution of the last coupling
// iteration
auto& x0 = NumLib::GlobalVectorProvider::provider.getVector(x);
MathLib::LinAlg::copy(x, x0);
// append a solution vector of suitable size
solutions_of_last_cpl_iteration.emplace_back(&x0);
}
}
}
void ControlNode::setInitialConditions(double const t0)
{
for (auto& process_data : per_process_data)
{
auto& pcs = process_data->process;
auto const process_id = process_data->process_id;
auto& time_disc = *process_data->time_disc;
auto& ode_sys = *process_data->tdisc_ode_sys;
// append a solution vector of suitable size
process_solutions.emplace_back(
&NumLib::GlobalVectorProvider::provider.getVector(
ode_sys.getMatrixSpecifications(process_id)));
process_solutions_prev.emplace_back(
&NumLib::GlobalVectorProvider::provider.getVector(
ode_sys.getMatrixSpecifications(process_id)));
auto& x = *process_solutions[process_id];
auto& x_prev = *process_solutions_prev[process_id];
pcs.setInitialConditions(process_id, t0, x);
MathLib::LinAlg::finalizeAssembly(x);
time_disc.setInitialState(t0); // push IC
MathLib::LinAlg::copy(x, x_prev); // pushState
}
}
NumLib::NonlinearSolverStatus
ControlNode::solveCoupledEquationSystemsByStaggeredScheme(
const double t, const double dt, const std::size_t timestep_id,
std::unique_ptr<Output>& output)
{
// Coupling iteration
if (coupling_max_iterations != 0)
{
// Set the flag of the first iteration be true.
for (auto& conv_crit : coupling_conv_crit)
{
conv_crit->preFirstIteration();
}
}
auto resetCouplingConvergenceCriteria = [&]() {
for (auto& conv_crit : coupling_conv_crit)
{
conv_crit->reset();
}
};
NumLib::NonlinearSolverStatus nonlinear_solver_status{false, -1};
bool coupling_iteration_converged = true;
for (int global_coupling_iteration = 0;
global_coupling_iteration < coupling_max_iterations;
global_coupling_iteration++, resetCouplingConvergenceCriteria())
{
// TODO(wenqing): use process name
coupling_iteration_converged = true;
int const last_process_id = per_process_data.size() - 1;
for (auto& process_data : per_process_data)
{
auto const process_id = process_data->process_id;
BaseLib::RunTime time_timestep_process;
time_timestep_process.start();
// The following setting of coupled_solutions can be removed only if
// the CoupledSolutionsForStaggeredScheme and related functions are
// removed totally from the computation of the secondary variable
// and from post-time functions.
CoupledSolutionsForStaggeredScheme coupled_solutions(
process_solutions);
process_data->process.setCoupledSolutionsForStaggeredScheme(
&coupled_solutions);
nonlinear_solver_status = solveOneTimeStepOneProcess(
process_solutions, process_solutions_prev, timestep_id, t, dt,
*process_data, *output);
process_data->nonlinear_solver_status = nonlinear_solver_status;
INFO(
"[time] Solving process #{:d} took {:g} s in time step #{:d} "
" coupling iteration #{:d}",
process_id, time_timestep_process.elapsed(), timestep_id,
global_coupling_iteration);
if (!nonlinear_solver_status.error_norms_met)
{
WARN(
"The nonlinear solver failed in time step #{:d} at t = "
"{:g} s for process #{:d}.",
timestep_id, t, process_id);
last_step_rejected = true;
return nonlinear_solver_status;
}
// Check the convergence of the coupling iteration
auto& x = *process_solutions[process_id];
auto& x_old = *solutions_of_last_cpl_iteration[process_id];
if (global_coupling_iteration > 0)
{
MathLib::LinAlg::axpy(x_old, -1.0, x); // save dx to x_old
if (process_id == last_process_id)
{
INFO(
"------- Checking convergence criterion for coupled "
"solution -------");
coupling_conv_crit[process_id]->checkDeltaX(x_old, x);
coupling_iteration_converged =
coupling_iteration_converged &&
coupling_conv_crit[process_id]->isSatisfied();
}
}
MathLib::LinAlg::copy(x, x_old);
} // end of for (auto& process_data : _per_process_data)
if (coupling_iteration_converged && global_coupling_iteration > 0)
{
break;
}
if (!nonlinear_solver_status.error_norms_met)
{
return nonlinear_solver_status;
}
}
if (!coupling_iteration_converged)
{
WARN(
"The coupling iterations reaches its maximum number in time step"
"#{:d} at t = {:g} s",
timestep_id, t);
}
// if (_chemical_solver_interface)
// {
// // Sequential non-iterative approach applied here to perform water
// // chemistry calculation followed by resolving component transport
// // process.
// // TODO: move into a global loop to consider both mass balance
// over
// // space and localized chemical equilibrium between solutes.
// BaseLib::RunTime time_phreeqc;
// time_phreeqc.start();
// auto& pcs = _per_process_data[0]->process;
// _chemical_solver_interface->doWaterChemistryCalculation(
// pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions),
// dt);
// nonlinear_solver_status = solveOneTimeStepOneProcess(
// _process_solutions, _process_solutions_prev,
// timestep_id, t, dt,
// *_per_process_data[1], *_output);
// pcs.extrapolateIntegrationPointValuesToNodes(
// t,
// _chemical_solver_interface->getIntPtProcessSolutions(),
// _process_solutions);
// INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
// }
return nonlinear_solver_status;
}
static constexpr std::string_view timestepper_cannot_reduce_dt =
"Time stepper cannot reduce the time step size further.";
NumLib::NonlinearSolverStatus ControlNode::solveUncoupledEquationSystems(
const double t, const double dt, const std::size_t timestep_id,
std::unique_ptr<Output>& output)
{
NumLib::NonlinearSolverStatus nonlinear_solver_status;
for (auto& process_data : per_process_data)
{
auto const process_id = process_data->process_id;
nonlinear_solver_status = solveMonolithicProcess(
t, dt, timestep_id, *process_data, process_solutions,
process_solutions_prev, *output);
process_data->nonlinear_solver_status = nonlinear_solver_status;
if (!nonlinear_solver_status.error_norms_met)
{
ERR("The nonlinear solver failed in time step #{:d} at t = {:g} s "
"for process #{:d}.",
timestep_id, t, process_id);
if (!process_data->timestepper->canReduceTimestepSize())
{
// save unsuccessful solution
output->doOutputAlways(process_data->process, process_id,
timestep_id, t, process_solutions);
OGS_FATAL(timestepper_cannot_reduce_dt.data());
}
return nonlinear_solver_status;
}
}
return nonlinear_solver_status;
}
} // namespace ProcessLib
/**
* \file
*
* \copyright
* Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#pragma once
#include <memory>
#include <vector>
#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
#include "ProcessLib/ProcessData.h"
namespace NumLib
{
class ConvergenceCriterion;
}
namespace ProcessLib
{
class Output;
struct ControlNode
{
explicit ControlNode(
std::vector<std::unique_ptr<ProcessData>>&& per_process_data_,
int const coupling_max_iterations_,
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>&&
coupling_conv_criteria_)
: per_process_data(std::move(per_process_data_)),
coupling_max_iterations(coupling_max_iterations_),
coupling_conv_crit(std::move(coupling_conv_criteria_))
{
}
void calculateNonEquilibriumInitialResiduum();
NumLib::NonlinearSolverStatus solveEquationSystems(
const double t,
const double dt,
const std::size_t timestep_id,
std::unique_ptr<Output>& output);
bool isMonolithicProcess()
{
// All _per_process_data share the first process.
return per_process_data[0]->process.isMonolithicSchemeUsed();
}
/**
* This function fills the vector of solutions of coupled processes of
* processes, _solutions_of_coupled_processes, and initializes the vector
* of solutions of the previous coupling iteration,
* _solutions_of_last_cpl_iteration.
*/
void setCoupledSolutions();
/**
* This function fills the vector of solutions of coupled processes of
* processes, _solutions_of_coupled_processes, and initializes the vector
* of solutions of the previous coupling iteration,
* _solutions_of_last_cpl_iteration.
*/
void setInitialConditions(double const t0);
/**
* \brief Member to solver coupled systems of equations by the staggered
* scheme.
*
* @param t Current time
* @param dt Time step size
* @param timestep_id Index of the time step
* @return true: if all nonlinear solvers convergence.
* false: if any of nonlinear solvers divergences.
*/
NumLib::NonlinearSolverStatus solveCoupledEquationSystemsByStaggeredScheme(
const double t,
const double dt,
const std::size_t timestep_id,
std::unique_ptr<Output>& output);
/**
* \brief Member to solver non coupled systems of equations, which can be
* a single system of equations, or several systems of equations
* without any dependency among the different systems.
*
* @param t Current time
* @param dt Time step size
* @param timestep_id Index of the time step
* @return true: if all nonlinear solvers convergence.
* false: if any of nonlinear solvers divergences.
*/
NumLib::NonlinearSolverStatus solveUncoupledEquationSystems(
const double t,
const double dt,
const std::size_t timestep_id,
std::unique_ptr<Output>& output);
std::vector<std::unique_ptr<ProcessData>> per_process_data;
const int coupling_max_iterations;
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>
coupling_conv_crit;
bool last_step_rejected = false;
std::vector<GlobalVector*> process_solutions;
std::vector<GlobalVector*> process_solutions_prev;
/// Solutions of the previous coupling iteration for the convergence
/// criteria of the coupling iteration.
std::vector<GlobalVector*> solutions_of_last_cpl_iteration;
};
} // namespace ProcessLib
......@@ -14,9 +14,9 @@
#ifdef USE_PETSC
#include "NumLib/ODESolver/PETScNonlinearSolver.h"
#endif // USE_PETSC
#include "NumLib/TimeStepping/CreateTimeStepper.h"
#include "ControlNode.h"
#include "CreateProcessData.h"
#include "NumLib/TimeStepping/CreateTimeStepper.h"
namespace ProcessLib
{
......@@ -64,18 +64,14 @@ static std::unique_ptr<ProcessData> makeProcessData(
OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
}
std::vector<std::vector<std::unique_ptr<ProcessData>>> createPerProcessData(
std::vector<ControlNode> createPerProcessData(
BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<Process>> const& processes,
std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>> const&
nonlinear_solvers)
{
std::vector<std::vector<std::unique_ptr<ProcessData>>> process_data;
std::vector<ControlNode> control_nodes;
int process_id = 0;
auto start_time = std::numeric_limits<double>::quiet_NaN();
auto end_time = std::numeric_limits<double>::quiet_NaN();
//! \ogs_file_param{prj__time_loop__processes__process}
for (auto pcs_config : config.getConfigSubtreeList("process"))
{
......@@ -88,6 +84,36 @@ std::vector<std::vector<std::unique_ptr<ProcessData>>> createPerProcessData(
},
"A process with the given name has not been defined.");
auto const& coupling_config
//! \ogs_file_param{prj__time_loop__global_process_coupling}
= pcs_config.getConfigSubtreeOptional("process_coupling");
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>
coupling_conv_criteria;
int max_coupling_iterations = 1;
if (coupling_config)
{
max_coupling_iterations
//! \ogs_file_param{prj__time_loop__global_process_coupling__max_iter}
= coupling_config->getConfigParameter<int>("max_iter");
auto const& coupling_convergence_criteria_config =
//! \ogs_file_param{prj__time_loop__global_process_coupling__convergence_criteria}
coupling_config->getConfigSubtree("convergence_criteria");
for (
auto coupling_convergence_criterion_config :
//! \ogs_file_param{prj__time_loop__global_process_coupling__convergence_criteria__convergence_criterion}
coupling_convergence_criteria_config.getConfigSubtreeList(
"convergence_criterion"))
{
coupling_conv_criteria.push_back(
NumLib::createConvergenceCriterion(
coupling_convergence_criterion_config));
}
}
int process_id = 0;
std::vector<std::unique_ptr<ProcessData>> per_process_data;
for (auto pv_config :
pcs_config.getConfigSubtreeList("process_variable"))
......@@ -147,9 +173,11 @@ std::vector<std::vector<std::unique_ptr<ProcessData>>> createPerProcessData(
"staggered scheme.");
}
process_data.push_back(std::move(per_process_data));
control_nodes.emplace_back(std::move(per_process_data),
max_coupling_iterations,
std::move(coupling_conv_criteria));
}
return process_data;
return control_nodes;
}
} // namespace ProcessLib
......@@ -14,7 +14,9 @@
namespace ProcessLib
{
std::vector<std::vector<std::unique_ptr<ProcessData>>> createPerProcessData(
struct ControlNode;
std::vector<ControlNode> createPerProcessData(
BaseLib::ConfigTree const& config,