Skip to content

GitLab

Commit 46a92d1a authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

worked on 1D problems

parent 6a05716d
Pipeline #75631 failed with stages
in 71 minutes and 14 seconds
Hide whitespace changes
Inline Side-by-side
mlmc/conf/m++.conf
View file @ 46a92d1a
#loadconf = mlmc_transport.conf
loadconf = mlmc_elliptic.conf
\ No newline at end of file
loadconf = mlmc_transport.conf
#loadconf = mlmc_elliptic.conf
\ No newline at end of file
mlmc/conf/mlmc_elliptic.conf
View file @ 46a92d1a
Model = LagrangeFEM
#Model = LagrangeFEM
#Model = MixedFEM
#Model = HybridFEM
Model = HybridFEM
#Model = DGFEM
degree = 1
......@@ -8,8 +8,9 @@ degree = 1
#sign = 1
#penalty = 20
#Problem = StochasticLaplace1D
Problem = StochasticLaplace2D
Problem = StochasticLaplace1D
#Problem = DeterministicLaplace1D
#Problem = StochasticLaplace2D
functional = L2
#functional = Energy
......@@ -41,7 +42,7 @@ MLMCVerbose = 2
PDEVerbose = 0
GeneratorVerbose = 0
GeneratorPlotting = 0
GeneratorPlotting = 1
MCPlotting = 0
MLMCPlotting = 0
......
mlmc/conf/mlmc_transport.conf
View file @ 46a92d1a
......@@ -3,14 +3,14 @@ Model = DGTransport
Overlap = dG1
flux_alpha = 1 #Upwind: 1, Central: 0
#Problem = StochasticPollution1D
Problem = StochasticPollution1D
#Problem = DeterministicPollution1D
#Problem = StochasticPollution2D
#Problem = DeterministicPollution2D
#Problem = StochasticPollutionCosHat2D
#Problem = DeterministicPollutionCosHat2D
#Problem = DeterministicPollutionCosHatSquare500
Problem = StochasticPollutionMollifiedBar2D
#Problem = StochasticPollutionMollifiedBar2D
#Problem = DeterministicPollutionMollifiedBar2D
scaling = 8
TimeSeries = uniform
......@@ -36,8 +36,8 @@ StochasticField = LogNormal
#StochasticField = Gaussian
#Experiment = ConvergenceTest
#Experiment = MLMCExperiment
Experiment = MLMCOverEpsilon
Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon
initLevels = 4, 5, 6
initSampleAmount = 16, 4, 2
......@@ -47,8 +47,9 @@ epsilon = 0.01
mcOnly = false
epsilonList = 0.01, 0.005, 0.003, 0.001
maxLevel = 9
maxLevel = 7
plevel = 3
degree = 2
MainVerbose = 1
MCVerbose = 2
......
mlmc/src/MLMCMain.C
View file @ 46a92d1a
......@@ -12,7 +12,6 @@ void MLMCMain::initialize() {
assemble = getAssemble();
plots = new MultilevelPlotter(*meshes);
// auto mgContainer = new MatrixGraphContainer(*meshes);
mlmc = new MultilevelMonteCarlo(initLevels, initSampleAmount, meshes, *plots,
stochasticField, assemble);
......@@ -63,9 +62,11 @@ StochasticField *MLMCMain::getStochasticField() {
}
if (problemName.find("Deterministic") != string::npos) {
if (problemName.find("Laplace") != string::npos)
return nullptr;
return new StochasticField(*meshes,
"DeterministicPermeabilityGenerator");
if (problemName.find("Pollution") != string::npos)
return new StochasticField(*meshes, "DeterministicHybridFluxGenerator");
return new StochasticField(*meshes,
"DeterministicHybridFluxGenerator");
}
}
......@@ -89,10 +90,11 @@ Assemble *MLMCMain::getAssemble() {
}
void MLMCMain::Run() {
printInfo();
this->printInfo();
mlmc->PrintInfo();
meshes->PrintInfo();
// assemble->PrintInfo();
// Stochastic field info
if (experimentName == "ConvergenceTest") {
headConvergenceTest();
......
mlmc/src/problem/StochasticProblemFactory.hpp
View file @ 46a92d1a
......@@ -19,7 +19,7 @@ enum EStochasticProblem {
DeterministicPollution1D,
DeterministicPollution2D,
DeterministicPollutionCosHat2D,
DeterministicPollutionCosHatSquare500
DeterministicPollutionMollifiedBar2D
};
class StochasticProblemFactory {
......@@ -41,9 +41,9 @@ public:
return (T *) new class StochasticPollution2D();
case StochasticPollutionCosHat2D:
case DeterministicPollutionCosHat2D:
case DeterministicPollutionCosHatSquare500:
return (T *) new class StochasticPollutionCosHat2D();
case StochasticPollutionMollifiedBar2D:
case DeterministicPollutionMollifiedBar2D:
return (T *) new class StochasticPollutionMollifiedBar2D();
}
}
......@@ -68,8 +68,8 @@ private:
{"DeterministicPollution2D", DeterministicPollution2D},
{"StochasticPollutionCosHat2D", StochasticPollutionCosHat2D},
{"DeterministicPollutionCosHat2D", DeterministicPollutionCosHat2D},
{"DeterministicPollutionCosHatSquare500",
DeterministicPollutionCosHatSquare500},
{"DeterministicPollutionMollifiedBar2D",
DeterministicPollutionMollifiedBar2D},
{"StochasticPollutionMollifiedBar2D", StochasticPollutionMollifiedBar2D}
};
}
......
mlmc/src/problem/StochasticTransportProblem.hpp
View file @ 46a92d1a
......@@ -6,8 +6,15 @@
class StochasticTransportProblem : public StochasticProblem {
protected:
Discretization *disc = nullptr;
HybridEllipticAssemble *assemble = nullptr;
public:
StochasticTransportProblem() = default;
~StochasticTransportProblem() override {
delete disc;
delete assemble;
}
VectorField b = Point(0.0, -1.0, 0.0);
......@@ -24,21 +31,38 @@ public:
virtual VectorField B(const Point &x) const { return b; }
virtual VectorField cellB(const cell &c, const Point &x) const {
return B(x);
return assemble->EvaluateCellFlux(*fieldSample, c);
}
virtual double faceB(const cell &c, int f, const VectorField &N,
const Point &x) const {
return B(x) * N;
virtual double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const {
return assemble->EvaluateNormalFlux(*fieldSample, c, f);
}
};
class StochasticPollution1D : public StochasticTransportProblem {
private:
double c = 1.0;
public:
StochasticPollution1D() = default;
StochasticPollution1D() {
disc = new Discretization("RT0_P0", 1, 1);
assemble = new HybridEllipticAssemble(disc, nullptr);
}
double ut(double t, const Point &x) const override {
return 0.0;
if (abs(c * t - x[0] + 0.5) > 0.06251) return 0.0;
return 1.0;
}
VectorField cellB(const cell &c, const Point &x) const override {
return VectorField(1.0, 0.0);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return VectorField(1.0, 0.0) * N;
}
string Name() const override { return "Stochastic Pollution 1D"; }
......@@ -46,7 +70,10 @@ public:
class StochasticPollution2D : public StochasticTransportProblem {
public:
StochasticPollution2D() = default;
StochasticPollution2D() {
disc = new Discretization("RT0_P0", 1, 2);
assemble = new HybridEllipticAssemble(disc, nullptr);
}
double ut(double t, const Point &x) const override {
double d = 1.0 / 16.0;
......@@ -55,28 +82,14 @@ public:
else return 0.0;
}
VectorField cellB(const cell &c, const Point &x) const override {
VectorField localFlux;
for (int d = 0; d < 3; d++)
localFlux[d] = (*fieldSample)(c(), d);
return localFlux;
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return N * VectorField(0.0, -1.0);
}
string Name() const override { return "Stochastic Pollution 2D"; }
};
class StochasticPollutionCosHat2D : public StochasticTransportProblem {
private:
Discretization *disc;
HybridEllipticAssemble *assemble;
double aa = 1.00;
double amplitude = 1.00;
double cc = 6.0;
public:
StochasticPollutionCosHat2D() {
disc = new Discretization("RT0_P0", 1, 2);
......@@ -87,69 +100,52 @@ public:
Point midPoint = Point(0.5, 0.8);
double r = dist(midPoint, x);
if (r < 1 / cc)
return aa * pow(cos(cc * Pi * r) + 1.0, 2.0);
return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
return 0.0;
}
VectorField cellB(const cell &c, const Point &x) const override {
return assemble->EvaluateCellFlux (*fieldSample, c);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return assemble->EvaluateNormalFlux (*fieldSample, c, f);
}
string Name() const override { return "Stochastic Pollution 2D"; }
};
class StochasticPollutionMollifiedBar2D : public StochasticTransportProblem{
class StochasticPollutionMollifiedBar2D : public StochasticTransportProblem {
private:
double ee = 0.12;
double xl = 0.6;
double yl = 0.05;
Point midPoint = Point(0.5,0.8);
double sh = 1/(xl*yl + 2*(xl+yl)*exp(1)*ee*0.2219969080840397 + exp(1)*2*Pi*ee*ee*0.0742477533879610);
Point tl = Point(midPoint[0]-0.5*xl,midPoint[1]+0.5*yl);
Point tr = Point(midPoint[0]+0.5*xl,midPoint[1]+0.5*yl);
Point br = Point(midPoint[0]+0.5*xl,midPoint[1]-0.5*yl);
Point bl = Point(midPoint[0]-0.5*xl,midPoint[1]-0.5*yl);
Point midPoint = Point(0.5, 0.8);
double sh = 1 / (xl * yl + 2 * (xl + yl) * exp(1) * ee * 0.2219969080840397
+ exp(1) * 2 * Pi * ee * ee * 0.0742477533879610);
Point tl = Point(midPoint[0] - 0.5 * xl, midPoint[1] + 0.5 * yl);
Point tr = Point(midPoint[0] + 0.5 * xl, midPoint[1] + 0.5 * yl);
Point br = Point(midPoint[0] + 0.5 * xl, midPoint[1] - 0.5 * yl);
Point bl = Point(midPoint[0] - 0.5 * xl, midPoint[1] - 0.5 * yl);
double phit(const Point &x, const Point &y) const {
if (dist(x,y)<ee)
return sh*exp(1.0) *exp(-1 / (1 - pow(dist(x,y) / ee,2)));
if (dist(x, y) < ee)
return sh * exp(1.0) * exp(-1 / (1 - pow(dist(x, y) / ee, 2)));
return 0.0;
}
Discretization *disc;
HybridEllipticAssemble *assemble;
public:
StochasticPollutionMollifiedBar2D(){
disc = new Discretization("RT0_P0",1,2);
StochasticPollutionMollifiedBar2D() {
disc = new Discretization("RT0_P0", 1, 2);
assemble = new HybridEllipticAssemble(disc, nullptr);
}
double ut(double t, const Point &x) const override {
if (dist(midPoint[0],x[0])<= 0.5*xl && dist(midPoint[1],x[1])<=0.5*yl) {
if (dist(midPoint[0], x[0]) <= 0.5 * xl && dist(midPoint[1], x[1]) <= 0.5 * yl) {
return sh;
}else if(dist(midPoint[0],x[0])<= 0.5*xl) {
} else if (dist(midPoint[0], x[0]) <= 0.5 * xl) {
Point hp1 = Point(x[0], midPoint[1] - 0.5 * yl);
Point hp2 = Point(x[0], midPoint[1] + 0.5 * yl);
return max(this->phit(x, hp1),this->phit(x, hp2));
}else if(dist(midPoint[1],x[1])<=0.5*yl){
Point hp1 = Point(midPoint[0]-0.5*xl,x[1]);
Point hp2 = Point(midPoint[0]+0.5*xl,x[1]);
return max(this->phit(x,hp1),this->phit(x,hp2));
}else
return max(max(this->phit(x,tl) ,this->phit(x,tr)),max(this->phit(x,br) ,this->phit(x,bl)));
}
VectorField cellB(const cell &c, const Point &x) const override {
return assemble->EvaluateCellFlux (*fieldSample, c);
}
double faceB(const cell &c, int f,
const VectorField &N, const Point &x) const override {
return assemble->EvaluateNormalFlux (*fieldSample, c, f);
return max(this->phit(x, hp1), this->phit(x, hp2));
} else if (dist(midPoint[1], x[1]) <= 0.5 * yl) {
Point hp1 = Point(midPoint[0] - 0.5 * xl, x[1]);
Point hp2 = Point(midPoint[0] + 0.5 * xl, x[1]);
return max(this->phit(x, hp1), this->phit(x, hp2));
} else
return max(max(this->phit(x, tl), this->phit(x, tr)),
max(this->phit(x, br), this->phit(x, bl)));
}
string Name() const override { return "Stochastic Pollution 2D"; }
......
mlmc/src/stochastics/StochasticField.h
View file @ 46a92d1a
......@@ -43,6 +43,12 @@ private:
if (meshes.dim() == 2)
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 (meshes.dim() == 1)
return new HybridFluxGenerator1D(l, meshes);
......
python/test_experiments.py
View file @ 46a92d1a
......@@ -61,11 +61,18 @@ class TestEllipticConvergence(TestExperiments):
self.assertEqual(len(os.listdir(mpp.PROJECT_PY_DATA_DIR)), 6)
class TestTransportTransport(TestExperiments):
def test_run(self):
epsilons = [0.05, 0.03, 0.01, 0.005]
for epsilon in epsilons:
kwargs = {'epsilon': epsilon}
class TestTransportBasics(TestExperiments):
def test_problems(self):
problems = [#'StochasticPollution1D',
#'DeterministicPollution1D',
'StochasticPollution2D',
'DeterministicPollution2D',
'StochasticPollutionCosHat2D',
'DeterministicPollutionCosHat2D',
'StochasticPollutionMollifiedBar2D',
'DeterministicPollutionMollifiedBar2D']
for problem in problems:
kwargs = {'Problem': problem}
self.assertEqual(mpp.run(4, config='mlmc_transport', kwargs=kwargs), 0)
def test_file_creation(self):
......
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

Imprint