Skip to content
Snippets Groups Projects
Commit 7171ec02 authored by Norihiro Watanabe's avatar Norihiro Watanabe
Browse files

changes after removing Mat/Vec templates from NumLib/ODESolver

parent 4e7010d7
No related branches found
No related tags found
No related merge requests found
......@@ -166,7 +166,7 @@ void ProjectData::buildProcesses()
_nonlinear_solvers, nl_slv_name,
"A nonlinear solver with the given name has not been defined.");
auto time_disc = NumLib::createTimeDiscretization<GlobalVector>(
auto time_disc = NumLib::createTimeDiscretization(
//! \ogs_file_param{process__time_discretization}
pc.getConfigSubtree("time_discretization"));
......
......@@ -40,7 +40,6 @@ class Mesh;
namespace NumLib
{
template <typename Matrix, typename Vector>
class NonlinearSolverBase;
}
......@@ -199,8 +198,7 @@ private:
std::unique_ptr<MathLib::LinearSolver<GlobalMatrix, GlobalVector>>>
_linear_solvers;
using NonlinearSolver =
NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>;
using NonlinearSolver = NumLib::NonlinearSolverBase;
std::map<std::string, std::unique_ptr<NonlinearSolver>> _nonlinear_solvers;
std::map<std::string,
std::unique_ptr<MathLib::PiecewiseLinearInterpolation>>
......
......@@ -39,14 +39,13 @@ public:
private:
//! An abstract nonlinear solver
using AbstractNLSolver =
NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>;
using AbstractNLSolver = NumLib::NonlinearSolverBase;
//! An abstract equations system
using EquationSystem = NumLib::EquationSystem<GlobalVector>;
using EquationSystem = NumLib::EquationSystem;
//! An abstract process
using Process = ProcessLib::Process;
//! An abstract time discretization
using TimeDisc = NumLib::TimeDiscretization<GlobalVector>;
using TimeDisc = NumLib::TimeDiscretization;
std::vector<GlobalVector*> _process_solutions;
std::unique_ptr<NumLib::ITimeStepAlgorithm> _timestepper;
......@@ -54,15 +53,12 @@ private:
struct SingleProcessData
{
template <NumLib::ODESystemTag ODETag, NumLib::NonlinearSolverTag NLTag>
SingleProcessData(NumLib::NonlinearSolver<GlobalMatrix, GlobalVector,
NLTag>& nonlinear_solver,
SingleProcessData(NumLib::NonlinearSolver<NLTag>& nonlinear_solver,
TimeDisc& time_disc,
NumLib::ODESystem<GlobalMatrix, GlobalVector, ODETag,
NLTag>& ode_sys)
NumLib::ODESystem<ODETag, NLTag>& ode_sys)
: nonlinear_solver_tag(NLTag),
nonlinear_solver(nonlinear_solver),
tdisc_ode_sys(new NumLib::TimeDiscretizedODESystem<
GlobalMatrix, GlobalVector, ODETag, NLTag>(
tdisc_ode_sys(new NumLib::TimeDiscretizedODESystem<ODETag, NLTag>(
ode_sys, time_disc)),
mat_strg(
dynamic_cast<NumLib::InternalMatrixStorage&>(*tdisc_ode_sys))
......@@ -104,17 +100,14 @@ private:
template <NumLib::ODESystemTag ODETag>
SingleProcessData makeSingleProcessData(
AbstractNLSolver& nonlinear_solver,
NumLib::ODESystem<GlobalMatrix, GlobalVector, ODETag,
NumLib::NonlinearSolverTag::Picard>& ode_sys,
NumLib::ODESystem<ODETag, NumLib::NonlinearSolverTag::Picard>& ode_sys,
TimeDisc& time_disc)
{
using Tag = NumLib::NonlinearSolverTag;
// A concrete Picard solver
using NonlinearSolverPicard =
NumLib::NonlinearSolver<GlobalMatrix, GlobalVector, Tag::Picard>;
using NonlinearSolverPicard = NumLib::NonlinearSolver<Tag::Picard>;
// A concrete Newton solver
using NonlinearSolverNewton =
NumLib::NonlinearSolver<GlobalMatrix, GlobalVector, Tag::Newton>;
using NonlinearSolverNewton = NumLib::NonlinearSolver<Tag::Newton>;
if (auto* nonlinear_solver_picard =
dynamic_cast<NonlinearSolverPicard*>(&nonlinear_solver))
......@@ -131,8 +124,7 @@ private:
{
// The Newton-Raphson method needs a Newton-ready ODE.
using ODENewton = NumLib::ODESystem<GlobalMatrix, GlobalVector,
ODETag, Tag::Newton>;
using ODENewton = NumLib::ODESystem<ODETag, Tag::Newton>;
if (auto* ode_newton = dynamic_cast<ODENewton*>(&ode_sys))
{
return SingleProcessData{*nonlinear_solver_newton, time_disc,
......@@ -172,10 +164,8 @@ private:
static void setEquationSystem(AbstractNLSolver& nonlinear_solver,
EquationSystem& eq_sys)
{
using Solver =
NumLib::NonlinearSolver<GlobalMatrix, GlobalVector, NLTag>;
using EqSys =
NumLib::NonlinearSystem<GlobalMatrix, GlobalVector, NLTag>;
using Solver = NumLib::NonlinearSolver<NLTag>;
using EqSys = NumLib::NonlinearSystem<NLTag>;
assert(dynamic_cast<Solver*>(&nonlinear_solver) != nullptr);
assert(dynamic_cast<EqSys*>(&eq_sys) != nullptr);
......@@ -210,8 +200,7 @@ private:
static void applyKnownSolutions(EquationSystem const& eq_sys,
GlobalVector& x)
{
using EqSys =
NumLib::NonlinearSystem<GlobalMatrix, GlobalVector, NLTag>;
using EqSys = NumLib::NonlinearSystem<NLTag>;
assert(dynamic_cast<EqSys const*>(&eq_sys) != nullptr);
auto& eq_sys_ = static_cast<EqSys const&>(eq_sys);
......
......@@ -103,7 +103,7 @@ public:
*
* \param timeDisc the time discretization scheme to be used.
*/
MatrixTranslatorGeneral(TimeDiscretization<Vector> const& timeDisc)
MatrixTranslatorGeneral(TimeDiscretization const& timeDisc)
: _time_disc(timeDisc)
{
}
......@@ -200,7 +200,7 @@ public:
*
* \param timeDisc the time discretization scheme to be used.
*/
MatrixTranslatorForwardEuler(ForwardEuler<Vector> const& timeDisc)
MatrixTranslatorForwardEuler(ForwardEuler const& timeDisc)
: _fwd_euler(timeDisc)
{
}
......
......@@ -17,37 +17,37 @@
namespace NumLib
{
template <typename Vector>
std::unique_ptr<TimeDiscretization<Vector>> createTimeDiscretization(
std::unique_ptr<TimeDiscretization> createTimeDiscretization(
BaseLib::ConfigTree const& config)
{
using T = std::unique_ptr<TimeDiscretization<Vector>>;
using T = std::unique_ptr<TimeDiscretization>;
//! \ogs_file_param{process__time_discretization__type}
auto const type = config.getConfigParameter<std::string>("type");
if (type == "BackwardEuler")
{
using ConcreteTD = BackwardEuler<Vector>;
using ConcreteTD = BackwardEuler;
return T(new ConcreteTD);
}
else if (type == "ForwardEuler")
{
using ConcreteTD = ForwardEuler<Vector>;
using ConcreteTD = ForwardEuler;
return T(new ConcreteTD);
}
else if (type == "CrankNicolson")
{
//! \ogs_file_param{process__time_discretization__CrankNicolson__theta}
auto const theta = config.getConfigParameter<double>("theta");
using ConcreteTD = CrankNicolson<Vector>;
using ConcreteTD = CrankNicolson;
return T(new ConcreteTD(theta));
}
else if (type == "BackwardDifferentiationFormula")
{
//! \ogs_file_param{process__time_discretization__BackwardDifferentiationFormula__order}
auto const order = config.getConfigParameter<unsigned>("order");
using ConcreteTD = BackwardDifferentiationFormula<Vector>;
using ConcreteTD = BackwardDifferentiationFormula;
return T(new ConcreteTD(order));
}
else
......
......@@ -31,7 +31,7 @@ GroundwaterFlowProcess::GroundwaterFlowProcess(
std::move(process_output)),
_process_data(std::move(process_data))
{
if (dynamic_cast<NumLib::ForwardEuler<GlobalVector>*>(
if (dynamic_cast<NumLib::ForwardEuler*>(
&Base::getTimeDiscretization()) != nullptr)
{
OGS_FATAL(
......
......@@ -25,17 +25,15 @@ class Mesh;
namespace ProcessLib
{
class Process : public NumLib::ODESystem<
GlobalMatrix, GlobalVector,
// TODO: later on use a simpler ODE system
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear,
NumLib::NonlinearSolverTag::Newton>
class Process
: public NumLib::ODESystem< // TODO: later on use a simpler ODE system
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear,
NumLib::NonlinearSolverTag::Newton>
{
public:
using Index = GlobalMatrix::IndexType;
using NonlinearSolver =
NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>;
using TimeDiscretization = NumLib::TimeDiscretization<GlobalVector>;
using NonlinearSolver = NumLib::NonlinearSolverBase;
using TimeDiscretization = NumLib::TimeDiscretization;
Process(MeshLib::Mesh& mesh,
NonlinearSolver& nonlinear_solver,
......
......@@ -11,8 +11,8 @@
#include "ODEs.h"
using EDMatrix = Eigen::MatrixXd;
using EVector = Eigen::VectorXd;
//using EDMatrix = Eigen::MatrixXd;
//using EVector = Eigen::VectorXd;
using GMatrix = GlobalMatrix;
using GVector = GlobalVector;
......@@ -22,9 +22,9 @@ template<typename Matrix, typename Vector, NumLib::NonlinearSolverTag NLTag>
class TestOutput
{
public:
using TimeDisc = NumLib::TimeDiscretization<Vector>;
using TimeDisc = NumLib::TimeDiscretization;
using LinearSolver = MathLib::LinearSolver<Matrix, Vector>;
using NLSolver = NumLib::NonlinearSolver<Matrix, Vector, NLTag>;
using NLSolver = NumLib::NonlinearSolver<NLTag>;
explicit TestOutput(const char* name)
: _file_name_part(name)
......@@ -36,21 +36,22 @@ public:
run_test<ODE>(ode, timeDisc, 10); // by default make 10 timesteps
}
template<template<typename /*Matrix*/, typename /*Vector*/> class ODE>
void run_test(ODE<Matrix, Vector>& ode, TimeDisc& timeDisc, const unsigned num_timesteps)
template <class ODE>
void run_test(ODE& ode, TimeDisc& timeDisc, const unsigned num_timesteps)
{
using ODE_ = ODE<Matrix, Vector>;
using ODET = ODETraits<Matrix, Vector, ODE>;
using ODE_ = ODE;
using ODET = ODETraits<ODE>;
NumLib::TimeDiscretizedODESystem<Matrix, Vector, ODE_::ODETag, NLTag>
NumLib::TimeDiscretizedODESystem<ODE_::ODETag, NLTag>
ode_sys(ode, timeDisc);
auto linear_solver = MathLib::createLinearSolver<Matrix, Vector>(nullptr);
std::unique_ptr<NLSolver> nonlinear_solver(
new NLSolver(*linear_solver, _tol, _maxiter));
NumLib::TimeLoopSingleODE<Matrix, Vector, NLTag> loop(
ode_sys, std::move(linear_solver), std::move(nonlinear_solver));
NumLib::TimeLoopSingleODE<NLTag> loop(ode_sys, std::move(linear_solver),
std::move(nonlinear_solver));
const double t0 = ODET::t0;
const double t_end = ODET::t_end;
......@@ -97,10 +98,10 @@ private:
*_file << "\n";
}
template<template<typename /*Matrix*/, typename /*Vector*/> class Ode>
template <class Ode>
void loopCallback(const double t, Vector const& x)
{
write(t, x, ODETraits<Matrix, Vector, Ode>::solution(t));
write(t, x, ODETraits<Ode>::solution(t));
}
const std::string _file_name_part;
......@@ -113,7 +114,7 @@ private:
template<typename Matrix, typename Vector, typename TimeDisc, typename ODE,
NumLib::NonlinearSolverTag NLTag>
typename std::enable_if<std::is_same<TimeDisc, NumLib::BackwardEuler<Vector> >::value>::type
typename std::enable_if<std::is_same<TimeDisc, NumLib::BackwardEuler >::value>::type
run_test_case(const unsigned num_timesteps, const char* name)
{
ODE ode;
......@@ -125,7 +126,7 @@ run_test_case(const unsigned num_timesteps, const char* name)
template<typename Matrix, typename Vector, typename TimeDisc, typename ODE,
NumLib::NonlinearSolverTag NLTag>
typename std::enable_if<std::is_same<TimeDisc, NumLib::ForwardEuler<Vector> >::value>::type
typename std::enable_if<std::is_same<TimeDisc, NumLib::ForwardEuler >::value>::type
run_test_case(const unsigned num_timesteps, const char* name)
{
ODE ode;
......@@ -137,7 +138,7 @@ run_test_case(const unsigned num_timesteps, const char* name)
template<typename Matrix, typename Vector, typename TimeDisc, typename ODE,
NumLib::NonlinearSolverTag NLTag>
typename std::enable_if<std::is_same<TimeDisc, NumLib::CrankNicolson<Vector> >::value>::type
typename std::enable_if<std::is_same<TimeDisc, NumLib::CrankNicolson >::value>::type
run_test_case(const unsigned num_timesteps, const char* name)
{
ODE ode;
......@@ -150,7 +151,7 @@ run_test_case(const unsigned num_timesteps, const char* name)
template<typename Matrix, typename Vector, typename TimeDisc, typename ODE,
NumLib::NonlinearSolverTag NLTag>
typename std::enable_if<
std::is_same<TimeDisc, NumLib::BackwardDifferentiationFormula<Vector> >::value>::type
std::is_same<TimeDisc, NumLib::BackwardDifferentiationFormula >::value>::type
run_test_case(const unsigned num_timesteps, const char* name)
{
ODE ode;
......@@ -163,24 +164,19 @@ run_test_case(const unsigned num_timesteps, const char* name)
// This class is only here s.t. I don't have to put the members into
// the definition of the macro TCLITEM below.
template<typename Matrix_, typename Vector_,
template<typename /*Matrix*/, typename /*Vector*/> class ODE_,
template<typename /*Vector*/> class TimeDisc_,
NumLib::NonlinearSolverTag NLTag_>
template <typename Matrix_, typename Vector_, class ODE_, class TimeDisc_,
NumLib::NonlinearSolverTag NLTag_>
struct TestCaseBase
{
using Matrix = Matrix_;
using Vector = Vector_;
using ODE = ODE_<Matrix_, Vector_>;
using TimeDisc = TimeDisc_<Vector_>;
using ODE = ODE_;
using TimeDisc = TimeDisc_;
static const NumLib::NonlinearSolverTag NLTag = NLTag_;
};
template<typename Matrix_, typename Vector_,
template<typename /*Matrix*/, typename /*Vector*/> class ODE_,
template<typename /*Vector*/> class TimeDisc_,
NumLib::NonlinearSolverTag NLTag_>
template <typename Matrix_, typename Vector_, class ODE_, class TimeDisc_,
NumLib::NonlinearSolverTag NLTag_>
struct TestCase;
......@@ -189,28 +185,29 @@ struct TestCase;
// Put new test cases to that list
//
// /////////////////////////////////////
//#define TESTCASESLIST \
// /* Eigen dense matrix */ \
// TCLITEM(EDMatrix, EVector, ODE1, BackwardEuler, Newton) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE1, ForwardEuler, Newton) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE1, CrankNicolson, Newton) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE1, BackwardDifferentiationFormula, Newton) TCLSEP \
// \
// TCLITEM(EDMatrix, EVector, ODE1, BackwardEuler, Picard) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE1, ForwardEuler, Picard) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE1, CrankNicolson, Picard) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE1, BackwardDifferentiationFormula, Picard) TCLSEP \
// \
// TCLITEM(EDMatrix, EVector, ODE2, BackwardEuler, Newton) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE2, ForwardEuler, Newton) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE2, CrankNicolson, Newton) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE2, BackwardDifferentiationFormula, Newton) TCLSEP \
// \
// TCLITEM(EDMatrix, EVector, ODE2, BackwardEuler, Picard) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE2, ForwardEuler, Picard) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE2, CrankNicolson, Picard) TCLSEP \
// TCLITEM(EDMatrix, EVector, ODE2, BackwardDifferentiationFormula, Picard) TCLSEP \
//
#define TESTCASESLIST \
/* Eigen dense matrix */ \
TCLITEM(EDMatrix, EVector, ODE1, BackwardEuler, Newton) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE1, ForwardEuler, Newton) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE1, CrankNicolson, Newton) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE1, BackwardDifferentiationFormula, Newton) TCLSEP \
\
TCLITEM(EDMatrix, EVector, ODE1, BackwardEuler, Picard) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE1, ForwardEuler, Picard) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE1, CrankNicolson, Picard) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE1, BackwardDifferentiationFormula, Picard) TCLSEP \
\
TCLITEM(EDMatrix, EVector, ODE2, BackwardEuler, Newton) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE2, ForwardEuler, Newton) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE2, CrankNicolson, Newton) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE2, BackwardDifferentiationFormula, Newton) TCLSEP \
\
TCLITEM(EDMatrix, EVector, ODE2, BackwardEuler, Picard) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE2, ForwardEuler, Picard) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE2, CrankNicolson, Picard) TCLSEP \
TCLITEM(EDMatrix, EVector, ODE2, BackwardDifferentiationFormula, Picard) TCLSEP \
\
/* Global sparse matrix */ \
TCLITEM(GMatrix, GVector, ODE1, BackwardEuler, Newton) TCLSEP \
TCLITEM(GMatrix, GVector, ODE1, ForwardEuler, Newton) TCLSEP \
......@@ -308,17 +305,16 @@ TYPED_TEST(NumLibODEIntTyped, T1)
TEST(NumLibODEInt, ODE3)
{
const char* name = "dummy";
{
// {
// only make sure ODE3 compiles
run_test_case<EDMatrix, EVector, NumLib::BackwardEuler<EVector>,
ODE3<EDMatrix, EVector>,
NumLib::NonlinearSolverTag::Newton>(0u, name);
}
// // only make sure ODE3 compiles
// run_test_case<EDMatrix, EVector, NumLib::BackwardEuler,
// ODE3<EDMatrix, EVector>,
// NumLib::NonlinearSolverTag::Newton>(0u, name);
// }
{
run_test_case<GMatrix, GVector, NumLib::BackwardEuler<GVector>,
ODE3<GMatrix, GVector>,
run_test_case<GMatrix, GVector, NumLib::BackwardEuler, ODE3,
NumLib::NonlinearSolverTag::Newton>(0u, name);
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment