Commit d7f16f2b authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactored PDESolver

parent cdfd7ab1
Pipeline #158348 canceled with stages
in 3 minutes and 54 seconds
...@@ -33,20 +33,20 @@ void EllipticPDESolver::plotSolution(SampleSolution &solution) { ...@@ -33,20 +33,20 @@ void EllipticPDESolver::plotSolution(SampleSolution &solution) {
IDiscretization *pressureDisc; IDiscretization *pressureDisc;
if (typeid(*assemble).name() == typeid(HybridEllipticAssemble).name()) if (typeid(*assemble).name() == typeid(HybridEllipticAssemble).name())
pressureDisc = new LagrangeDiscretization(meshes, 0, 1); pressureDisc = new LagrangeDiscretization(GetDisc()->GetMeshes(), 0, 1);
else else
pressureDisc = new LagrangeDiscretization(meshes, 1, 1); pressureDisc = new LagrangeDiscretization(GetDisc()->GetMeshes(), 1, 1);
SampleSolution pressure(pressureDisc, solution.Id(), "Pressure"); SampleSolution pressure(pressureDisc, solution.Id(), "Pressure");
assemble->SetPressure(solution.U, pressure.U); assemble->SetPressure(solution.U, pressure.U);
mpp::plot(pressure.IdString()) << pressure.U << mpp::endp; mpp::plot(pressure.IdString()) << pressure.U << mpp::endp;
LagrangeDiscretization fluxDisc(meshes, 0, SpaceDimension); LagrangeDiscretization fluxDisc(GetDisc()->GetMeshes(), 0, SpaceDimension);
SampleSolution flux(&fluxDisc, solution.Id(), "Flux"); SampleSolution flux(&fluxDisc, solution.Id(), "Flux");
assemble->SetFlux(solution.U, flux.U); assemble->SetFlux(solution.U, flux.U);
mpp::plot(flux.IdString()) << flux.U << mpp::endp; mpp::plot(flux.IdString()) << flux.U << mpp::endp;
LagrangeDiscretization kappaDisc(meshes, 0, 1); LagrangeDiscretization kappaDisc(GetDisc()->GetMeshes(), 0, 1);
SampleSolution kappa(&kappaDisc, solution.Id(), "Kappa"); SampleSolution kappa(&kappaDisc, solution.Id(), "Kappa");
assemble->GetProblem()->Permeability(kappa.U); assemble->GetProblem()->Permeability(kappa.U);
mpp::plot(kappa.IdString()) << kappa.U << mpp::endp; mpp::plot(kappa.IdString()) << kappa.U << mpp::endp;
......
#ifndef ELLIPTICPDESOLVER_HPP #ifndef ELLIPTICPDESOLVER_HPP
#define ELLIPTICPDESOLVER_HPP #define ELLIPTICPDESOLVER_HPP
#include "LagrangeEllipticAssemble.hpp" #include "LagrangeEllipticAssemble.hpp"
#include "MixedEllipticAssemble.hpp" #include "MixedEllipticAssemble.hpp"
#include "HybridEllipticAssemble.hpp" #include "HybridEllipticAssemble.hpp"
...@@ -12,7 +11,11 @@ ...@@ -12,7 +11,11 @@
class EllipticPDESolver : public PDESolver { class EllipticPDESolver : public PDESolver {
private: private:
Newton *newton; std::unique_ptr<Newton> newton;
std::string quantity = "L2";
std::string costMeasure = "size";
IStochasticEllipticAssemble *assemble; IStochasticEllipticAssemble *assemble;
...@@ -26,21 +29,21 @@ protected: ...@@ -26,21 +29,21 @@ protected:
void plotSolution(SampleSolution &solution) override; void plotSolution(SampleSolution &solution) override;
public: public:
EllipticPDESolver(IStochasticEllipticAssemble *assemble, const Meshes &meshes, EllipticPDESolver(IStochasticEllipticAssemble *assemble,
const std::string &quantity = "L2", const string &quantity, const string &costMeasure) :
const std::string &costMeasure = "size") : PDESolver(), assemble(assemble), quantity(quantity), costMeasure(costMeasure),
PDESolver(meshes, quantity, costMeasure), newton(std::make_unique<Newton>()) {
newton(new Newton()), assemble(assemble) {
if (verbose > 0) { if (verbose > 0) {
PrintInfo(); mout.PrintInfo("EllipticPDESolver", verbose,
PrintInfoEntry("Quantity", quantity),
PrintInfoEntry("Cost measure", costMeasure));
newton->PrintInfo(); newton->PrintInfo();
assemble->PrintInfo(); assemble->PrintInfo();
} }
} }
~EllipticPDESolver() { ~EllipticPDESolver() override {
delete assemble; delete assemble;
delete newton;
} }
IAssemble *GetAssemble() const override { return assemble; } IAssemble *GetAssemble() const override { return assemble; }
......
...@@ -18,12 +18,6 @@ protected: ...@@ -18,12 +18,6 @@ protected:
int plotting = 0; int plotting = 0;
const Meshes &meshes;
std::string quantity;
std::string costMeasure;
virtual void run(SampleSolution &solution) = 0; virtual void run(SampleSolution &solution) = 0;
virtual void computeQ(SampleSolution &solution) = 0; virtual void computeQ(SampleSolution &solution) = 0;
...@@ -37,23 +31,13 @@ protected: ...@@ -37,23 +31,13 @@ protected:
} }
public: public:
PDESolver(const Meshes &meshes, PDESolver() {
const std::string &quantity,
const std::string &costMeasure) :
meshes(meshes), quantity(quantity), costMeasure(costMeasure) {
config.get("PDESolverVerbose", verbose); config.get("PDESolverVerbose", verbose);
config.get("PDESolverPlotting", plotting); config.get("PDESolverPlotting", plotting);
} }
virtual ~PDESolver() = default; virtual ~PDESolver() = default;
virtual void PrintInfo() const {
if (verbose > 0)
mout.PrintInfo(Name(), verbose,
PrintInfoEntry("Quantity", quantity),
PrintInfoEntry("Cost measure", costMeasure));
}
void Run(SampleSolution &solution) { void Run(SampleSolution &solution) {
mout.StartBlock("PDE Solver"); mout.StartBlock("PDE Solver");
vout(2) << solution.IdString() << endl; vout(2) << solution.IdString() << endl;
...@@ -95,9 +79,8 @@ protected: ...@@ -95,9 +79,8 @@ protected:
void plotSolution(SampleSolution &solution) override {} void plotSolution(SampleSolution &solution) override {}
public: public:
// Todo remove costMeasure and quantity explicit DummyPDESolver(IStochasticDummyAssemble *assemble) :
DummyPDESolver(IStochasticDummyAssemble *assemble, const Meshes &meshes) : PDESolver(), assemble(assemble) {}
PDESolver(meshes, "", ""), assemble(assemble) {}
IAssemble *GetAssemble() const override { return assemble; } IAssemble *GetAssemble() const override { return assemble; }
......
...@@ -6,42 +6,42 @@ ...@@ -6,42 +6,42 @@
PDESolver *PDESolverCreator::Create(const Meshes &meshes) { PDESolver *PDESolverCreator::Create(const Meshes &meshes) {
if (_model == "LagrangeElliptic") if (_model == "LagrangeElliptic")
return new EllipticPDESolver( return new EllipticPDESolver(
new LagrangeEllipticAssemble( new LagrangeEllipticAssemble(
new LagrangeDiscretization(meshes, _degree), new LagrangeDiscretization(meshes, _degree),
CreateStochasticEllipticProblem(_problem, meshes) CreateStochasticEllipticProblem(_problem, meshes)
), meshes, _quantity, _costMeasure ), _quantity, _costMeasure
); );
if (_model == "MixedElliptic") if (_model == "MixedElliptic")
return new EllipticPDESolver( return new EllipticPDESolver(
new MixedEllipticAssemble( new MixedEllipticAssemble(
new RTLagrangeDiscretization(meshes, 0, 0), new RTLagrangeDiscretization(meshes, 0, 0),
CreateStochasticEllipticProblem(_problem, meshes) CreateStochasticEllipticProblem(_problem, meshes)
), meshes, _quantity, _costMeasure ), _quantity, _costMeasure
); );
if (_model == "HybridElliptic") if (_model == "HybridElliptic")
return new EllipticPDESolver( return new EllipticPDESolver(
new HybridEllipticAssemble( new HybridEllipticAssemble(
new HybridRTLagrangeDiscretization(meshes, 0, 0), new HybridRTLagrangeDiscretization(meshes, 0, 0),
CreateStochasticEllipticProblem(_problem, meshes) CreateStochasticEllipticProblem(_problem, meshes)
), meshes, _quantity, _costMeasure ), _quantity, _costMeasure
); );
if (_model == "DGElliptic") if (_model == "DGElliptic")
return new EllipticPDESolver( return new EllipticPDESolver(
new DGEllipticAssemble( new DGEllipticAssemble(
new DGDiscretization(meshes, _degree), new DGDiscretization(meshes, _degree),
CreateStochasticEllipticProblem(_problem, meshes) CreateStochasticEllipticProblem(_problem, meshes)
), meshes, _quantity, _costMeasure ), _quantity, _costMeasure
); );
if (_model == "DGTransport") if (_model == "DGTransport")
return new TransportPDESolver( return new TransportPDESolver(
new DGLinearTransportAssemble( new DGLinearTransportAssemble(
new DGDiscretization(meshes, _degree), new DGDiscretization(meshes, _degree),
CreateStochasticTransportProblem(_problem, meshes) CreateStochasticTransportProblem(_problem, meshes)
), meshes, _quantity, _costMeasure ), _quantity, _costMeasure
); );
if (_model == "PGTransport") if (_model == "PGTransport")
...@@ -65,11 +65,10 @@ PDESolver *PDESolverCreator::Create(const Meshes &meshes) { ...@@ -65,11 +65,10 @@ PDESolver *PDESolverCreator::Create(const Meshes &meshes) {
if (_model == "DummyPDESolver") if (_model == "DummyPDESolver")
return new DummyPDESolver( return new DummyPDESolver(
new IStochasticDummyAssemble( new IStochasticDummyAssemble(
new LagrangeDiscretization(meshes, _degree), new LagrangeDiscretization(meshes, _degree),
CreateStochasticDummyProblem(_problem, meshes) CreateStochasticDummyProblem(_problem, meshes)
), meshes ));
);
Exit(_model + " not found") Exit(_model + " not found")
} }
......
...@@ -20,6 +20,9 @@ void TransportPDESolver::computeQ(SampleSolution &solution) { ...@@ -20,6 +20,9 @@ void TransportPDESolver::computeQ(SampleSolution &solution) {
} }
void TransportPDESolver::computeCost(SampleSolution &solution) { void TransportPDESolver::computeCost(SampleSolution &solution) {
if (costMeasure == "size") solution.C = solution.U.size();
// else if (costMeasure == "time") solution.Cost = solution.U.size(); // Todo
else Exit("Cost measure not implemented")
} }
void TransportPDESolver::plotSolution(SampleSolution &solution) { void TransportPDESolver::plotSolution(SampleSolution &solution) {
......
...@@ -8,7 +8,11 @@ ...@@ -8,7 +8,11 @@
class TransportPDESolver : public PDESolver { class TransportPDESolver : public PDESolver {
private: private:
TimeIntegrator *timeInt; std::string quantity = "L2";
std::string costMeasure = "size";
std::unique_ptr<TimeIntegrator> timeInt;
IStochasticLinearTransportAssemble *assemble; IStochasticLinearTransportAssemble *assemble;
...@@ -23,23 +27,22 @@ protected: ...@@ -23,23 +27,22 @@ protected:
public: public:
TransportPDESolver(IStochasticLinearTransportAssemble *assemble, TransportPDESolver(IStochasticLinearTransportAssemble *assemble,
const Meshes &meshes, const string &quantity, const string &costMeasure) :
const std::string &quantity = "L2", PDESolver(), assemble(assemble), quantity(quantity), costMeasure(costMeasure),
const std::string &costMeasure = "size") :
PDESolver(meshes, quantity, costMeasure), assemble(assemble),
timeInt(TimeIntegratorCreator(IMPLICIT_MIDPOINT). timeInt(TimeIntegratorCreator(IMPLICIT_MIDPOINT).
WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))). WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))).
WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))). WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))).
WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))). WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))).
WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))). WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))).
Create()) { CreateUnique()) {
if (verbose) if (verbose)
PrintInfo(); mout.PrintInfo("TransportPDESolver", verbose,
PrintInfoEntry("Quantity", quantity),
PrintInfoEntry("Cost measure", costMeasure));
} }
~TransportPDESolver() { ~TransportPDESolver() override {
delete assemble; delete assemble;
delete timeInt;
} }
IAssemble *GetAssemble() const override { IAssemble *GetAssemble() const override {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment