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

refactored HybridFluxGenerator, started with fixing transport

parent 2027608e
...@@ -21,7 +21,7 @@ public: ...@@ -21,7 +21,7 @@ public:
std::shared_ptr<StochasticTransportProblem> problem) std::shared_ptr<StochasticTransportProblem> problem)
: DGTAssemble(dynamic_cast<DGDiscretization *>(disc.get())), : DGTAssemble(dynamic_cast<DGDiscretization *>(disc.get())),
problem(move(problem)) { problem(move(problem)) {
ReadConfig(Settings, "flux_alpha", flux_alpha); config.get("flux_alpha", flux_alpha);
logger = IndentedLogger::GetInstance(); logger = IndentedLogger::GetInstance();
} }
......
...@@ -46,13 +46,12 @@ void MonteCarloTransport::method() { ...@@ -46,13 +46,12 @@ void MonteCarloTransport::method() {
void MonteCarloTransport::solvePDE(int l, int m, bool fineLevel, void MonteCarloTransport::solvePDE(int l, int m, bool fineLevel,
double &valueQ, double &cost, Vector &solution) { double &valueQ, double &cost, Vector &solution) {
logger->StartMethod("Start Solving PDE", verbose); logger->StartMethod("Start Solving PDE", verbose);
logger->InnerMsg("l: " + to_string(l) + logger->InnerMsg("l: " + to_string(l) + " m: " + to_string(m), verbose);
" m: " + to_string(m), verbose);
Preconditioner *pc = GetPC("PointBlockGaussSeidel"); //TODO RMV after usage Preconditioner *pc = GetPC("PointBlockGaussSeidel");
Solver solver(pc, "GMRES"); Solver solver(pc, "GMRES");
Solver solver2(pc, "GMRES"); // TODO RMV
DGTimeIntegrator timeIntegrator(solver, solver, solver, solver2); DGTimeIntegrator timeIntegrator(solver, solver, solver, solver);
TimeSeries timeSeries = getTimeSeries(fineLevel); TimeSeries timeSeries = getTimeSeries(fineLevel);
assemble->resetTAssemble(); assemble->resetTAssemble();
......
...@@ -29,6 +29,9 @@ private: ...@@ -29,6 +29,9 @@ private:
public: public:
DGTransportAssemble *assemble; DGTransportAssemble *assemble;
unique_ptr<Preconditioner> pc;
unique_ptr<Solver> solver;
MatrixGraphs cellMatrixGraphs; MatrixGraphs cellMatrixGraphs;
MatrixGraphs faceMatrixGraphs; MatrixGraphs faceMatrixGraphs;
CellMatrixGraphs solMatrixGraphs; CellMatrixGraphs solMatrixGraphs;
...@@ -38,11 +41,8 @@ public: ...@@ -38,11 +41,8 @@ public:
Vector fineSolution; Vector fineSolution;
Vector coarseSolution; Vector coarseSolution;
MonteCarloTransport(int l, MonteCarloTransport(int l, int dM, bool baseLevel,
int dM, Meshes *meshes, StochasticField *stochFields,
bool baseLevel,
Meshes *meshes,
StochasticField *stochFields,
DGTransportAssemble *assemble) : DGTransportAssemble *assemble) :
MonteCarlo(l, dM, baseLevel, meshes, stochFields), MonteCarlo(l, dM, baseLevel, meshes, stochFields),
assemble(assemble), assemble(assemble),
...@@ -52,7 +52,9 @@ public: ...@@ -52,7 +52,9 @@ public:
fineInput(Vector(faceMatrixGraphs[l - pLevel])), fineInput(Vector(faceMatrixGraphs[l - pLevel])),
coarseInput(Vector(faceMatrixGraphs[l - pLevel - 1])), coarseInput(Vector(faceMatrixGraphs[l - pLevel - 1])),
fineSolution(Vector(solMatrixGraphs[l - pLevel])), fineSolution(Vector(solMatrixGraphs[l - pLevel])),
coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])) { coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])),
pc(unique_ptr<Preconditioner>(GetPC("PointBlockGaussSeidel"))),
solver(make_unique<Solver>(pc.get(), "GMRES")) {
config.get("startTime", startTime); config.get("startTime", startTime);
config.get("endTime", endTime); config.get("endTime", endTime);
......
...@@ -5,8 +5,9 @@ ...@@ -5,8 +5,9 @@
class StochasticProblem { class StochasticProblem {
public: protected:
Vector *fieldSample = nullptr; Vector *fieldSample = nullptr;
public:
StochasticProblem() = default; StochasticProblem() = default;
......
...@@ -5,13 +5,9 @@ void HybridFluxGenerator::generateFineSample(Vector &fineField) { ...@@ -5,13 +5,9 @@ void HybridFluxGenerator::generateFineSample(Vector &fineField) {
initialize(true); initialize(true);
permeabilityGenerator->SetSampleDir(sampleDir); permeabilityGenerator->SetSampleDir(sampleDir);
permeabilityGenerator->GetFineSample(*finePerm); permeabilityGenerator->GetFineSample(*finePerm);
assemble->problem->SetSample(finePerm); assemble->problem->SetSample(finePerm.get());
Preconditioner *pc = GetPC("SuperLU"); newton(*assemble, *u);
Solver solver(pc, "GMRES");
Newton newton(solver);
(newton)(*assemble, *u);
assemble->setFlux(*u, *flux); assemble->setFlux(*u, *flux);
assemble->SetNormalFlux(*u, fineField); assemble->SetNormalFlux(*u, fineField);
...@@ -24,13 +20,9 @@ void HybridFluxGenerator::generateCoarseSample(const Vector &fineField, ...@@ -24,13 +20,9 @@ void HybridFluxGenerator::generateCoarseSample(const Vector &fineField,
Vector &coarseField) { Vector &coarseField) {
initialize(false); initialize(false);
permeabilityGenerator->GetCoarseSample(*finePerm, *coarsePerm); permeabilityGenerator->GetCoarseSample(*finePerm, *coarsePerm);
assemble->problem->SetSample(coarsePerm); assemble->problem->SetSample(coarsePerm.get());
Preconditioner *pc = GetPC("SuperLU");
Solver solver(pc, "GMRES");
Newton newton(solver);
(newton)(*assemble, *u); newton(*assemble, *u);
assemble->setFlux(*u, *flux); assemble->setFlux(*u, *flux);
assemble->SetNormalFlux(*u, coarseField); assemble->SetNormalFlux(*u, coarseField);
...@@ -41,11 +33,11 @@ void HybridFluxGenerator::generateCoarseSample(const Vector &fineField, ...@@ -41,11 +33,11 @@ void HybridFluxGenerator::generateCoarseSample(const Vector &fineField,
void HybridFluxGenerator::initialize(bool fineField) { void HybridFluxGenerator::initialize(bool fineField) {
if (fineField) { if (fineField) {
flux = new Vector(cellMatrixGraphs[l - meshes.pLevel()]); flux = make_unique<Vector>(cellMatrixGraphs[l - meshes.pLevel()]);
u = new Vector(faceMatrixGraphs[l - meshes.pLevel()]); u = make_unique<Vector>(faceMatrixGraphs[l - meshes.pLevel()]);
} else { } else {
flux = new Vector(cellMatrixGraphs[l - meshes.pLevel() - 1]); flux = make_unique<Vector>(cellMatrixGraphs[l - meshes.pLevel() - 1]);
u = new Vector(faceMatrixGraphs[l - meshes.pLevel() - 1]); u = make_unique<Vector>(faceMatrixGraphs[l - meshes.pLevel() - 1]);
} }
*flux = 0.0, *u = 0.0; *flux = 0.0, *u = 0.0;
} }
...@@ -12,45 +12,39 @@ private: ...@@ -12,45 +12,39 @@ private:
MatrixGraphs cellMatrixGraphs; MatrixGraphs cellMatrixGraphs;
MatrixGraphs faceMatrixGraphs; MatrixGraphs faceMatrixGraphs;
Vector *finePerm = nullptr; unique_ptr<Vector> finePerm;
Vector *coarsePerm = nullptr; unique_ptr<Vector> coarsePerm;
Vector *flux = nullptr; unique_ptr<Vector> flux;
Vector *u = nullptr; unique_ptr<Vector> u;
void initialize(bool fineField); void initialize(bool fineField);
protected: protected:
Discretization *disc = nullptr; shared_ptr<IDiscretizationT<>> disc;
StochasticEllipticProblem *problem = nullptr; shared_ptr<StochasticEllipticProblem> problem;
HybridEllipticAssemble *assemble = nullptr; unique_ptr<HybridEllipticAssemble> assemble;
CirculantEmbedding *permeabilityGenerator = nullptr; unique_ptr<CirculantEmbedding> permeabilityGenerator;
unique_ptr<Preconditioner> pc;
Solver solver;
Newton newton;
void generateFineSample(Vector &fineField) override; void generateFineSample(Vector &fineField) override;
void generateCoarseSample(const Vector &fineField, Vector &coarseField) override; void generateCoarseSample(const Vector &fineField, Vector &coarseField) override;
public: public:
explicit HybridFluxGenerator(int l, Meshes &meshes) explicit HybridFluxGenerator(int l, Meshes &meshes) :
: SampleGenerator(l, meshes), SampleGenerator(l, meshes),
cellMatrixGraphs(MatrixGraphs(meshes, dof("cell", 3))), cellMatrixGraphs(MatrixGraphs(meshes, dof("cell", 3))),
faceMatrixGraphs(MatrixGraphs(meshes, dof("face", 1))) { faceMatrixGraphs(MatrixGraphs(meshes, dof("face", 1))),
finePerm = new Vector(cellMatrixGraphs[l - meshes.pLevel()]); pc(unique_ptr<Preconditioner>(GetPC("SuperLU"))),
if (l != meshes.pLevel()) solver(Solver(pc.get(), "GMRES")),
coarsePerm = new Vector(cellMatrixGraphs[l - meshes.pLevel() - 1]); newton(Newton(solver)) {
else
coarsePerm = nullptr;
}
~HybridFluxGenerator() {
delete disc;
delete problem;
delete assemble;
delete permeabilityGenerator;
delete finePerm; finePerm = make_unique<Vector>(cellMatrixGraphs[l - meshes.pLevel()]);
delete coarsePerm; if (l != meshes.pLevel())
delete flux; coarsePerm = make_unique<Vector>(cellMatrixGraphs[l - meshes.pLevel() - 1]);
delete u;
} }
}; };
...@@ -58,11 +52,11 @@ class HybridFluxGenerator1D : public HybridFluxGenerator { ...@@ -58,11 +52,11 @@ class HybridFluxGenerator1D : public HybridFluxGenerator {
public: public:
explicit HybridFluxGenerator1D(int l, Meshes &meshes) : explicit HybridFluxGenerator1D(int l, Meshes &meshes) :
HybridFluxGenerator(l, meshes) { HybridFluxGenerator(l, meshes) {
disc = new Discretization("RT0_P0", 1, meshes.dim());
// problem = StochasticProblemFactory::CreateFrom<StochasticEllipticProblem>( disc = make_shared<DiscretizationT<>>(meshes, "RT0_P0");
// EStochasticProblem::StochasticLaplace1D); problem = make_shared<StochasticLaplace1D>();
// assemble = new HybridEllipticAssemble(disc, problem); assemble = make_unique<HybridEllipticAssemble>(disc, problem);
permeabilityGenerator = new CirculantEmbedding1D(l, meshes); permeabilityGenerator = make_unique<CirculantEmbedding1D>(l, meshes);
} }
}; };
...@@ -70,35 +64,11 @@ class HybridFluxGenerator2D : public HybridFluxGenerator { ...@@ -70,35 +64,11 @@ class HybridFluxGenerator2D : public HybridFluxGenerator {
public: public:
explicit HybridFluxGenerator2D(int l, Meshes &meshes) explicit HybridFluxGenerator2D(int l, Meshes &meshes)
: HybridFluxGenerator(l, meshes) { : HybridFluxGenerator(l, meshes) {
disc = new Discretization("RT0_P0", 1, meshes.dim());
// problem = StochasticProblemFactory::CreateFrom<StochasticEllipticProblem>(
// EStochasticProblem::StochasticLaplace2D);
// assemble = new HybridEllipticAssemble(disc, problem);
permeabilityGenerator = new CirculantEmbedding2D(l, meshes);
}
};
class DeterministicHybridFluxGenerator1D : public HybridFluxGenerator { disc = make_shared<DiscretizationT<>>(meshes, "RT0_P0");
public: problem = make_shared<StochasticLaplace2D>();
explicit DeterministicHybridFluxGenerator1D(int l, Meshes &meshes) assemble = make_unique<HybridEllipticAssemble>(disc, problem);
: HybridFluxGenerator(l, meshes) { permeabilityGenerator = make_unique<CirculantEmbedding2D>(l, meshes);
disc = new Discretization("RT0_P0", 1, meshes.dim());
// problem = StochasticProblemFactory::CreateFrom<StochasticEllipticProblem>(
// EStochasticProblem::StochasticLaplace1D);
// assemble = new HybridEllipticAssemble(disc, problem);
permeabilityGenerator = new DeterministicPermeabilityGenerator(l, meshes);
}
};
class DeterministicHybridFluxGenerator2D : public HybridFluxGenerator {
public:
explicit DeterministicHybridFluxGenerator2D(int l, Meshes &meshes)
: HybridFluxGenerator(l, meshes) {
disc = new Discretization("RT0_P0", 1, meshes.dim());
// problem = StochasticProblemFactory::CreateFrom<StochasticEllipticProblem>(
// EStochasticProblem::StochasticLaplace2D);
// assemble = new HybridEllipticAssemble(disc, problem);
permeabilityGenerator = new DeterministicPermeabilityGenerator(l, meshes);
} }
}; };
......
...@@ -35,24 +35,13 @@ private: ...@@ -35,24 +35,13 @@ private:
if (meshes.dim() == 2) if (meshes.dim() == 2)
return new CirculantEmbedding2D(l, meshes); return new CirculantEmbedding2D(l, meshes);
} }
if (generatorName == "DeterministicPermeabilityGenerator") {
// if (meshes.dim() == 1)
// return new DeterministicPermeabilityGenerator1D(l, meshes);
// if (meshes.dim() == 2)
// return new DeterministicPermeabilityGenerator2D(l, meshes);
}
if (generatorName == "HybridFluxGenerator") { if (generatorName == "HybridFluxGenerator") {
if (meshes.dim() == 1) if (meshes.dim() == 1)
return new HybridFluxGenerator1D(l, meshes); return new HybridFluxGenerator1D(l, meshes);
if (meshes.dim() == 2) if (meshes.dim() == 2)
return new HybridFluxGenerator2D(l, meshes); return new HybridFluxGenerator2D(l, meshes);
} }
if (generatorName == "DeterministicHybridFluxGenerator") {
if (meshes.dim() == 1)
return new DeterministicHybridFluxGenerator1D(l, meshes);
if (meshes.dim() == 2)
return new DeterministicHybridFluxGenerator2D(l, meshes);
}
Exit("Generator does not exist") Exit("Generator does not exist")
} }
}; };
......
...@@ -72,15 +72,15 @@ INSTANTIATE_TEST_CASE_P(TestMainProgram, TestMainProgram, Values( ...@@ -72,15 +72,15 @@ INSTANTIATE_TEST_CASE_P(TestMainProgram, TestMainProgram, Values(
std::map<std::string, std::string>{ std::map<std::string, std::string>{
{"GeneratorPlotting", "1"}, {"GeneratorPlotting", "1"},
{"MCPlotting", "1"}, {"MCPlotting", "1"},
{"Problem", "StochasticLaplace2D"}}}, {"Problem", "StochasticLaplace2D"}}}
// Tests with default transport config // Tests with default transport config
ConfigMapsForTest{defaultTransportConfigMap, // ConfigMapsForTest{defaultTransportConfigMap,
std::map<std::string, std::string>{}}, // std::map<std::string, std::string>{}}
//
ConfigMapsForTest{defaultTransportConfigMap, // ConfigMapsForTest{defaultTransportConfigMap,
std::map<std::string, std::string>{ // std::map<std::string, std::string>{
{"Problem", "StochasticPollution2D"}}} // {"Problem", "StochasticPollution2D"}}}
)); ));
TEST_P(TestMainProgram, TestInitialize) { TEST_P(TestMainProgram, TestInitialize) {
......
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