Commit 4bc3650f authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Merge branch '37-remove-meshes-from-pdesolver' into 'feature'

Resolve "Remove Meshes From PDESolver"

Closes #37

See merge request !50
parents cdfd7ab1 6046d4ac
Pipeline #159204 failed with stages
in 20 minutes and 8 seconds
......@@ -2,6 +2,7 @@ add_library(PDESOLVERS STATIC
PDESolverCreator.cpp
EllipticPDESolver.cpp
TransportPDESolver.cpp
ParabolicPDESolver.cpp
assembling/elliptic/LagrangeEllipticAssemble.cpp
assembling/elliptic/MixedEllipticAssemble.cpp
assembling/elliptic/HybridEllipticAssemble.cpp
......
......@@ -33,20 +33,20 @@ void EllipticPDESolver::plotSolution(SampleSolution &solution) {
IDiscretization *pressureDisc;
if (typeid(*assemble).name() == typeid(HybridEllipticAssemble).name())
pressureDisc = new LagrangeDiscretization(meshes, 0, 1);
pressureDisc = new LagrangeDiscretization(GetDisc()->GetMeshes(), 0, 1);
else
pressureDisc = new LagrangeDiscretization(meshes, 1, 1);
pressureDisc = new LagrangeDiscretization(GetDisc()->GetMeshes(), 1, 1);
SampleSolution pressure(pressureDisc, solution.Id(), "Pressure");
assemble->SetPressure(solution.U, pressure.U);
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");
assemble->SetFlux(solution.U, flux.U);
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");
assemble->GetProblem()->Permeability(kappa.U);
mpp::plot(kappa.IdString()) << kappa.U << mpp::endp;
......
#ifndef ELLIPTICPDESOLVER_HPP
#define ELLIPTICPDESOLVER_HPP
#include "LagrangeEllipticAssemble.hpp"
#include "MixedEllipticAssemble.hpp"
#include "HybridEllipticAssemble.hpp"
......@@ -12,7 +11,11 @@
class EllipticPDESolver : public PDESolver {
private:
Newton *newton;
std::unique_ptr<Newton> newton;
std::string quantity = "L2";
std::string costMeasure = "size";
IStochasticEllipticAssemble *assemble;
......@@ -26,21 +29,23 @@ protected:
void plotSolution(SampleSolution &solution) override;
public:
EllipticPDESolver(IStochasticEllipticAssemble *assemble, const Meshes &meshes,
const std::string &quantity = "L2",
const std::string &costMeasure = "size") :
PDESolver(meshes, quantity, costMeasure),
newton(new Newton()), assemble(assemble) {
EllipticPDESolver(IStochasticEllipticAssemble *assemble,
const string &quantity, const string &costMeasure) :
PDESolver(), assemble(assemble), quantity(quantity), costMeasure(costMeasure) {
newton = std::make_unique<Newton>();
if (verbose > 0) {
PrintInfo();
mout.PrintInfo("EllipticPDESolver", verbose,
PrintInfoEntry("Quantity", quantity),
PrintInfoEntry("Cost measure", costMeasure));
newton->PrintInfo();
assemble->PrintInfo();
}
}
~EllipticPDESolver() {
~EllipticPDESolver() override {
delete assemble;
delete newton;
}
IAssemble *GetAssemble() const override { return assemble; }
......
......@@ -18,12 +18,6 @@ protected:
int plotting = 0;
const Meshes &meshes;
std::string quantity;
std::string costMeasure;
virtual void run(SampleSolution &solution) = 0;
virtual void computeQ(SampleSolution &solution) = 0;
......@@ -37,23 +31,13 @@ protected:
}
public:
PDESolver(const Meshes &meshes,
const std::string &quantity,
const std::string &costMeasure) :
meshes(meshes), quantity(quantity), costMeasure(costMeasure) {
PDESolver() {
config.get("PDESolverVerbose", verbose);
config.get("PDESolverPlotting", plotting);
}
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) {
mout.StartBlock("PDE Solver");
vout(2) << solution.IdString() << endl;
......@@ -95,9 +79,8 @@ protected:
void plotSolution(SampleSolution &solution) override {}
public:
// Todo remove costMeasure and quantity
DummyPDESolver(IStochasticDummyAssemble *assemble, const Meshes &meshes) :
PDESolver(meshes, "", ""), assemble(assemble) {}
explicit DummyPDESolver(IStochasticDummyAssemble *assemble) :
PDESolver(), assemble(assemble) {}
IAssemble *GetAssemble() const override { return assemble; }
......
#include "PDESolverCreator.hpp"
#include "EllipticPDESolver.hpp"
#include "TransportPDESolver.hpp"
#include "ParabolicPDESolver.hpp"
PDESolver *PDESolverCreator::Create(const Meshes &meshes) {
if (_model == "LagrangeElliptic")
return new EllipticPDESolver(
new LagrangeEllipticAssemble(
new LagrangeDiscretization(meshes, _degree),
CreateStochasticEllipticProblem(_problem, meshes)
), meshes, _quantity, _costMeasure
new LagrangeEllipticAssemble(
new LagrangeDiscretization(meshes, _degree),
CreateStochasticEllipticProblem(_problem, meshes)
), _quantity, _costMeasure
);
if (_model == "MixedElliptic")
return new EllipticPDESolver(
new MixedEllipticAssemble(
new RTLagrangeDiscretization(meshes, 0, 0),
CreateStochasticEllipticProblem(_problem, meshes)
), meshes, _quantity, _costMeasure
new MixedEllipticAssemble(
new RTLagrangeDiscretization(meshes, 0, 0),
CreateStochasticEllipticProblem(_problem, meshes)
), _quantity, _costMeasure
);
if (_model == "HybridElliptic")
return new EllipticPDESolver(
new HybridEllipticAssemble(
new HybridRTLagrangeDiscretization(meshes, 0, 0),
CreateStochasticEllipticProblem(_problem, meshes)
), meshes, _quantity, _costMeasure
new HybridEllipticAssemble(
new HybridRTLagrangeDiscretization(meshes, 0, 0),
CreateStochasticEllipticProblem(_problem, meshes)
), _quantity, _costMeasure
);
if (_model == "DGElliptic")
return new EllipticPDESolver(
new DGEllipticAssemble(
new DGDiscretization(meshes, _degree),
CreateStochasticEllipticProblem(_problem, meshes)
), meshes, _quantity, _costMeasure
new DGEllipticAssemble(
new DGDiscretization(meshes, _degree),
CreateStochasticEllipticProblem(_problem, meshes)
), _quantity, _costMeasure
);
if (_model == "DGTransport")
return new TransportPDESolver(
new DGLinearTransportAssemble(
new DGDiscretization(meshes, _degree),
CreateStochasticTransportProblem(_problem, meshes)
), meshes, _quantity, _costMeasure
new DGLinearTransportAssemble(
new DGDiscretization(meshes, _degree),
CreateStochasticTransportProblem(_problem, meshes)
), _quantity, _costMeasure
);
if (_model == "PGTransport")
return nullptr; // Todo open research question
// if (_model == "DGReaction")
// return new ReactionPDESolver(
// new DGReactionAssemble(
// new DGDiscretization(meshes, _degree),
// CreateStochasticReactionProblem(_problem, meshes)
// ), meshes, _quantity, _costMeasure
// );
//
// if (_model == "PGReaction")
// return new ReactionPDESolver(
// new PGReactionAssemble(
// new DGDiscretization(meshes, _degree),
// CreateStochasticReactionProblem(_problem, meshes)
// ), meshes, _quantity, _costMeasure
// );
if (_model == "DGReaction")
return new ParabolicPDESolver(
new DGReactionAssemble(
new DGDiscretization(meshes, _degree),
CreateStochasticReactionProblem(_problem, meshes)
), _quantity, _costMeasure
);
if (_model == "PGReaction")
return new ParabolicPDESolver(
new PGReactionAssemble(
new DGDiscretization(meshes, _degree),
CreateStochasticReactionProblem(_problem, meshes)
), _quantity, _costMeasure
);
if (_model == "DummyPDESolver")
return new DummyPDESolver(
new IStochasticDummyAssemble(
new LagrangeDiscretization(meshes, _degree),
CreateStochasticDummyProblem(_problem, meshes)
), meshes
);
new IStochasticDummyAssemble(
new LagrangeDiscretization(meshes, _degree),
CreateStochasticDummyProblem(_problem, meshes)
));
Exit(_model + " not found")
}
......
//
// Created by niklas on 11.05.21.
//
#include "ParabolicPDESolver.hpp"
void ParabolicPDESolver::run(SampleSolution &solution) {
}
void ParabolicPDESolver::computeQ(SampleSolution &solution) {
}
void ParabolicPDESolver::computeCost(SampleSolution &solution) {
}
void ParabolicPDESolver::plotSolution(SampleSolution &solution) {
}
......@@ -4,7 +4,66 @@
#include "DGReactionAssemble.hpp"
#include "PGReactionAssemble.hpp"
class ParabolicPDESolver {
#include "PDESolver.hpp"
class ParabolicPDESolver : public PDESolver {
private:
std::string quantity = "L2";
std::string costMeasure = "size";
std::unique_ptr<NonLinearTimeIntegrator> timeInt;
IStochasticReactionAssemble *assemble;
protected:
void run(SampleSolution &solution) override;
void computeQ(SampleSolution &solution) override;
void computeCost(SampleSolution &solution) override;
void plotSolution(SampleSolution &solution) override;
public:
ParabolicPDESolver(IStochasticReactionAssemble *assemble,
const string &quantity, const string &costMeasure) :
PDESolver(), assemble(assemble), quantity(quantity), costMeasure(costMeasure) {
if (typeid(*assemble).name() == typeid(PGReactionAssemble).name())
timeInt = TimeIntegratorCreator(IMPLICIT_EULER).
WithNonLinearSolver(new Newton(std::make_unique<GMRES>(GetPC("SuperLU")))).
CreateUniqueNonLinearTimeIntegrator();
if (typeid(*assemble).name() == typeid(DGReactionAssemble).name())
timeInt = TimeIntegratorCreator(IMPLICIT_EULER).
WithNonLinearSolver(new Newton(std::make_unique<GMRES>(GetPC(
"PointBlockJacobi_dG")))).
CreateUniqueNonLinearTimeIntegrator();
if (verbose)
mout.PrintInfo("TransportPDESolver", verbose,
PrintInfoEntry("Quantity", quantity),
PrintInfoEntry("Cost measure", costMeasure));
}
~ParabolicPDESolver() override {
delete assemble;
}
IAssemble *GetAssemble() const override {
return nullptr;
// return assemble;
}
IDiscretization *GetDisc() const override { return nullptr; }
IStochasticProblem *GetProblem() const override { return nullptr; }
void DrawSample(const SampleID &id) override {}
std::string Name() const override { return "ParabolicPDESolver"; }
};
#endif //PARABOLICPDESOLVER_HPP
......@@ -14,12 +14,13 @@ void TransportPDESolver::computeQ(SampleSolution &solution) {
else if (quantity == "Error") solution.Q = assemble->Error(0, solution.U);
else if (quantity == "Inflow") solution.Q = assemble->InflowOutflow(solution.U).first;
else if (quantity == "Outflow") solution.Q = assemble->InflowOutflow(solution.U).second;
// else if (quantity == "L2Error") solution.Q = assemble->L2Error(solution.U);
// else if (quantity == "EnergyError") solution.Q = assemble->EnergyError(solution.U);
else Exit("Quantity of interest not implemented")
}
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) {
......
......@@ -8,7 +8,11 @@
class TransportPDESolver : public PDESolver {
private:
TimeIntegrator *timeInt;
std::string quantity = "L2";
std::string costMeasure = "size";
std::unique_ptr<TimeIntegrator> timeInt;
IStochasticLinearTransportAssemble *assemble;
......@@ -23,23 +27,24 @@ protected:
public:
TransportPDESolver(IStochasticLinearTransportAssemble *assemble,
const Meshes &meshes,
const std::string &quantity = "L2",
const std::string &costMeasure = "size") :
PDESolver(meshes, quantity, costMeasure), assemble(assemble),
timeInt(TimeIntegratorCreator(IMPLICIT_MIDPOINT).
const string &quantity, const string &costMeasure) :
PDESolver(), assemble(assemble), quantity(quantity), costMeasure(costMeasure) {
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"))).
Create()) {
CreateUnique();
if (verbose)
PrintInfo();
mout.PrintInfo("TransportPDESolver", verbose,
PrintInfoEntry("Quantity", quantity),
PrintInfoEntry("Cost measure", costMeasure));
}
~TransportPDESolver() {
~TransportPDESolver() override {
delete assemble;
delete timeInt;
}
IAssemble *GetAssemble() const override {
......
......@@ -16,7 +16,7 @@ public:
double flux_alpha = 1.0;
double diffusion = 0.0;
DGLinearTransportAssemble(IDiscretizationT<> *disc,
DGLinearTransportAssemble(IDiscretization *disc,
IStochasticTransportProblem *problem)
: IStochasticLinearTransportAssemble(problem),
disc(dynamic_cast<DGDiscretization *> (disc)) {
......
add_library(PROBLEMS STATIC
IStochasticProblem.cpp
StochasticEllipticProblem.cpp
StochasticReactionProblem.cpp
StochasticTransportProblem.cpp
)
target_link_libraries(PROBLEMS GENERATORS)
#include "StochasticEllipticProblem.hpp"
class StochasticLaplace1D : public IStochasticEllipticProblem {
CirculantEmbedding1D tensorGenerator;
public:
explicit StochasticLaplace1D(const Meshes &meshes) :
IStochasticEllipticProblem(meshes),
tensorGenerator(CirculantEmbedding1D(meshes)) {};
void DrawSample(const SampleID &id) override {
tensorGenerator.DrawSample(id);
}
Scalar Solution(const Point &x) const override {
return 1 - x[0];
}
VectorField Flux(const Point &x) const override {
return VectorField(1.0, 0.0);
}
Scalar Load(const Point &x) const override {
return 0.0;
}
Tensor Permeability(const cell &c) const override {
return this->tensorGenerator.EvalSample(c);
}
string Name() const override { return "StochasticLaplace1D"; }
};
class Laplace1D : public IStochasticEllipticProblem {
public:
explicit Laplace1D(const Meshes &meshes) :
IStochasticEllipticProblem(meshes) {};
void DrawSample(const SampleID &id) override {}
Scalar Solution(const Point &x) const override {
return 1 - x[0];
}
VectorField Flux(const Point &x) const override {
return VectorField(-1.0, 0.0);
}
Scalar Load(const Point &x) const override {
return 0.0;
}
Tensor Permeability(const cell &c) const override {
return One;
}
string Name() const override { return "Laplace1D"; }
};
class StochasticLaplace2D : public IStochasticEllipticProblem {
CirculantEmbedding2D tensorGenerator;
public:
explicit StochasticLaplace2D(const Meshes &meshes) :
IStochasticEllipticProblem(meshes),
tensorGenerator(CirculantEmbedding2D(meshes)) {}
void DrawSample(const SampleID &id) override {
tensorGenerator.DrawSample(id);
}
Scalar Solution(const Point &x) const override {
return -x[1];
}
VectorField Flux(const Point &x) const override {
return VectorField(0.0, -1.0);
}
Scalar Load(const Point &x) const override { return 0.0; }
Tensor Permeability(const cell &c) const override {
return this->tensorGenerator.EvalSample(c);
}
string Name() const override { return "StochasticLaplace2D"; }
};
class StochasticLaplace2DTest : public IStochasticEllipticProblem {
UniformDistributionRVector rVectorGenerator;
public:
explicit StochasticLaplace2DTest(const Meshes &meshes) :
IStochasticEllipticProblem(meshes),
rVectorGenerator(UniformDistributionRVector(meshes, 2, -1.0, 1.0)) {}
void DrawSample(const SampleID &id) override {
rVectorGenerator.DrawSample(id);
}
Scalar Solution(const Point &x) const override {
return -x[1];
}
VectorField Flux(const Point &x) const override {
return VectorField(0.0, -1.0);
}
Scalar Load(const Point &x) const override { return 0.0; }
Tensor Permeability(const cell &c) const override {
if (c()[1] > 0.125 && c()[1] < 0.375) {
if (c()[0] > 0.125 && c()[0] < 0.375)
return (1 + 0.9 * this->rVectorGenerator.EvalSample()[0]) * One;
if (c()[0] > 0.625 && c()[0] < 0.875)
return (1 + 0.9 * this->rVectorGenerator.EvalSample()[1]) * One;
}
return One;
}
string Name() const override { return "StochasticLaplace2DTest"; }
};
class Laplace2D : public IStochasticEllipticProblem {
public:
explicit Laplace2D(const Meshes &meshes) : IStochasticEllipticProblem(meshes) {}
void DrawSample(const SampleID &id) override {}
Scalar Solution(const Point &x) const override {
return -x[1];
}
VectorField Flux(const Point &x) const override {
return VectorField(0.0, -1.0);
}
Scalar Load(const Point &x) const override { return 0.0; }
Tensor Permeability(const cell &c) const override {
return One;
}
string Name() const override { return "Laplace2D"; }
};
IStochasticEllipticProblem *
CreateStochasticEllipticProblem(const string &problemName, const Meshes &meshes) {
if (problemName == "StochasticLaplace1D")
......
......@@ -23,147 +23,6 @@ public:
}
};
class StochasticLaplace1D : public IStochasticEllipticProblem {
CirculantEmbedding1D tensorGenerator;
public:
explicit StochasticLaplace1D(const Meshes &meshes) :
IStochasticEllipticProblem(meshes),
tensorGenerator(CirculantEmbedding1D(meshes)) {};
void DrawSample(const SampleID &id) override {
tensorGenerator.DrawSample(id);
}
Scalar Solution(const Point &x) const override {
return 1 - x[0];
}