Commit 19a9f39f authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

added 1D problems

parent 38acae70
#include "StochasticFields.h" #include "StochasticFields.h"
class StochasticFieldPair1D : public StochasticFieldPair {
int N;
CirculantEmbedding1D *circEmbedding1D;
vector<double> fieldf;
vector<double> fieldc;
public:
explicit StochasticFieldPair1D(int l) : StochasticFieldPair() {
// Todo rmv Hack!
N = (int) pow(2, l);
circEmbedding1D = new CirculantEmbedding1D(N);
}
~StochasticFieldPair1D() { delete circEmbedding1D; }
void generate() override {
if (PPM->master()) {
PPM->BroadcastInt(N);
fieldf = circEmbedding1D->getFieldf();
for (double &cell_value : fieldf) {
PPM->BroadcastDouble(cell_value);
}
fieldc = circEmbedding1D->generateInterpolField(fieldf);
for (double &cell_value : fieldc) {
PPM->BroadcastDouble(cell_value);
}
} else {
N = PPM->BroadcastInt();
fieldf.resize(N);
for (double &cell_value : fieldf)
cell_value = PPM->BroadcastDouble();
fieldc.resize(N / 2);
for (double &cell_value : fieldc)
cell_value = PPM->BroadcastDouble();
}
}
double getFineValue(const Point &z) const override {
int i = floor(z[0] * N);
return fieldf[i];
}
double getCoarseValue(const Point &z) const override {
int i = floor(z[0] * N / 2);
return fieldc[i];
}
pair<double, double> getFineMinMaxValues() const override {
double max = 0.0;
double min = infty;
for (int i = 0; i < N; i++) {
if (fieldf[i] > max) max = fieldf[i];
if (fieldf[i] < min) min = fieldf[i];
}
return {min, max};
}
pair<double, double> getCoarseMinMaxValues() const override {
double max = 0.0;
double min = infty;
for (int i = 0; i < N; i++) {
if (fieldc[i] > max) max = fieldc[i];
if (fieldc[i] < min) min = fieldc[i];
}
return {min, max};
}
};
class StochasticFieldPair2D : public StochasticFieldPair { class StochasticFieldPair2D : public StochasticFieldPair {
int N1; int N1;
int N2; int N2;
CirculantEmbedding2D *circEmbedding2D; CirculantEmbedding2D *circEmbedding2D;
std::vector<std::vector<double>> field_f; std::vector<std::vector<double>> fieldf;
std::vector<std::vector<double>> field_c; std::vector<std::vector<double>> fieldc;
public: public:
// explicit StochasticFieldPair2D(int l, const Mesh &coarseMesh, const Mesh &fineMesh) :
// StochasticFieldPair(coarseMesh, fineMesh) {
// // Todo rmv Hack!
// N1 = (int) pow(2, l);
// N2 = N1;
// circEmbedding2D = new CirculantEmbedding2D(N1, N2);
// }
explicit StochasticFieldPair2D(int l) : StochasticFieldPair() { explicit StochasticFieldPair2D(int l) : StochasticFieldPair() {
// Todo rmv Hack! // Todo rmv Hack!
N1 = (int) pow(2, l); N1 = (int) pow(2, l);
...@@ -28,28 +87,28 @@ public: ...@@ -28,28 +87,28 @@ public:
if (PPM->master()) { if (PPM->master()) {
PPM->BroadcastInt(N1); PPM->BroadcastInt(N1);
PPM->BroadcastInt(N2); PPM->BroadcastInt(N2);
field_f = circEmbedding2D->get_fieldf(); fieldf = circEmbedding2D->getFieldf();
for (int i = 0; i < N2; i++) for (int i = 0; i < N2; i++)
for (int j = 0; j < N1; j++) for (int j = 0; j < N1; j++)
PPM->BroadcastDouble(field_f[i][j]); PPM->BroadcastDouble(fieldf[i][j]);
field_c = circEmbedding2D->generate_interpol_field(field_f); fieldc = circEmbedding2D->generateInterpolField(fieldf);
for (int i = 0; i < N2 / 2; i++) for (int i = 0; i < N2 / 2; i++)
for (int j = 0; j < N1 / 2; j++) for (int j = 0; j < N1 / 2; j++)
PPM->BroadcastDouble(field_c[i][j]); PPM->BroadcastDouble(fieldc[i][j]);
} else { } else {
N1 = PPM->BroadcastInt(); N1 = PPM->BroadcastInt();
N2 = PPM->BroadcastInt(); N2 = PPM->BroadcastInt();
field_f.resize(N2); fieldf.resize(N2);
for (int i = 0; i < N2; i++) { for (int i = 0; i < N2; i++) {
field_f[i].resize(N1); fieldf[i].resize(N1);
for (int j = 0; j < N1; j++) for (int j = 0; j < N1; j++)
field_f[i][j] = PPM->BroadcastDouble(); fieldf[i][j] = PPM->BroadcastDouble();
} }
field_c.resize(N2 / 2); fieldc.resize(N2 / 2);
for (int i = 0; i < N2 / 2; i++) { for (int i = 0; i < N2 / 2; i++) {
field_c[i].resize(N1 / 2); fieldc[i].resize(N1 / 2);
for (int j = 0; j < N1 / 2; j++) for (int j = 0; j < N1 / 2; j++)
field_c[i][j] = PPM->BroadcastDouble(); fieldc[i][j] = PPM->BroadcastDouble();
} }
} }
} }
...@@ -57,13 +116,13 @@ public: ...@@ -57,13 +116,13 @@ public:
double getFineValue(const Point &z) const override { double getFineValue(const Point &z) const override {
int i = floor(z[0] * N2); int i = floor(z[0] * N2);
int j = floor(z[1] * N1); int j = floor(z[1] * N1);
return field_f[i][j]; return fieldf[i][j];
} }
double getCoarseValue(const Point &z) const override { double getCoarseValue(const Point &z) const override {
int i = floor(z[0] * N2 / 2); int i = floor(z[0] * N2 / 2);
int j = floor(z[1] * N1 / 2); int j = floor(z[1] * N1 / 2);
return field_c[i][j]; return fieldc[i][j];
} }
pair<double, double> getFineMinMaxValues() const override { pair<double, double> getFineMinMaxValues() const override {
...@@ -71,8 +130,8 @@ public: ...@@ -71,8 +130,8 @@ public:
double min = infty; double min = infty;
for (int i = 0; i < N2; i++) { for (int i = 0; i < N2; i++) {
for (int j = 0; j < N1; j++) { for (int j = 0; j < N1; j++) {
if (field_f[i][j] > max) max = field_f[i][j]; if (fieldf[i][j] > max) max = fieldf[i][j];
if (field_f[i][j] < min) min = field_f[i][j]; if (fieldf[i][j] < min) min = fieldf[i][j];
} }
} }
return {min, max}; return {min, max};
...@@ -83,21 +142,19 @@ public: ...@@ -83,21 +142,19 @@ public:
double min = infty; double min = infty;
for (int i = 0; i < N2 / 2; i++) { for (int i = 0; i < N2 / 2; i++) {
for (int j = 0; j < N1 / 2; j++) { for (int j = 0; j < N1 / 2; j++) {
if (field_c[i][j] > max) max = field_c[i][j]; if (fieldc[i][j] > max) max = fieldc[i][j];
if (field_c[i][j] < min) min = field_c[i][j]; if (fieldc[i][j] < min) min = fieldc[i][j];
} }
} }
return {min, max}; return {min, max};
} }
}; };
StochasticFieldPair *getStochasticFieldPair(int l, Meshes &meshes, StochasticFieldPair *getStochasticFieldPair(int l, int dim, string fieldName) {
string StochasticField) { if (dim == 1)
if (StochasticField.empty()) return new StochasticFieldPair1D(l);
ReadConfig(Settings, "StochasticField", StochasticField); else if (dim == 2)
if (StochasticField == "LogNormal2D")
// return new StochasticFieldPair2D(l, meshes[l - 1], meshes[l]);
return new StochasticFieldPair2D(l); return new StochasticFieldPair2D(l);
else Exit("This StochasticField is not implemented") else
return nullptr; return nullptr;
} }
\ No newline at end of file
...@@ -29,21 +29,19 @@ public: ...@@ -29,21 +29,19 @@ public:
}; };
}; };
StochasticFieldPair *getStochasticFieldPair(int l, Meshes &meshes, StochasticFieldPair *getStochasticFieldPair(int l, int dim, string fieldName);
string StochasticField = "");
class StochasticFields { class StochasticFields {
public: public:
const string &fieldName;
int currentLevel; int currentLevel;
bool fineField; bool fineField = true;
map<int, StochasticFieldPair *> stochFieldPairs; map<int, StochasticFieldPair *> stochFieldPairs;
explicit StochasticFields(Meshes &meshes) { explicit StochasticFields(const string &fieldName, Meshes &meshes) :
currentLevel = meshes.pLevel(); fieldName(fieldName), currentLevel(meshes.pLevel()) {
fineField = true;
for (int l = meshes.pLevel(); l <= meshes.Level(); l++) { for (int l = meshes.pLevel(); l <= meshes.Level(); l++) {
stochFieldPairs[l] = getStochasticFieldPair(l, meshes, stochFieldPairs[l] = getStochasticFieldPair(l, meshes.dim(), fieldName);
"LogNormal2D");
} }
} }
...@@ -63,8 +61,8 @@ public: ...@@ -63,8 +61,8 @@ public:
return stochFieldPairs[currentLevel]->getCoarseValue(z); return stochFieldPairs[currentLevel]->getCoarseValue(z);
}; };
virtual pair<double, double> getMinMaxValues() { virtual pair<double, double> getMinMaxValues(bool _fineField) {
if (fineField) if (_fineField)
return stochFieldPairs[currentLevel]->getFineMinMaxValues(); return stochFieldPairs[currentLevel]->getFineMinMaxValues();
else else
return stochFieldPairs[currentLevel]->getCoarseMinMaxValues(); return stochFieldPairs[currentLevel]->getCoarseMinMaxValues();
...@@ -74,10 +72,9 @@ public: ...@@ -74,10 +72,9 @@ public:
currentLevel = l; currentLevel = l;
} }
void setFineField(bool useFineField) { void setFineField(bool _fineField) {
useFineField = useFineField; fineField = _fineField;
} }
}; };
#endif //M_STOCHASTICFIELD_HPP #endif //M_STOCHASTICFIELD_HPP
#include "StochasticProblem.h" #include "StochasticProblem.h"
class StochasticProblem2D : public StochasticProblem { class StochasticLaplace1D : public StochasticProblem {
public: public:
explicit StochasticProblem2D(StochasticFields &stochFields) : explicit StochasticLaplace1D(StochasticFields &stochFields) :
StochasticProblem(stochFields) {}
Scalar Solution(const Point &x) const override { return -x[0]; }
VectorField Flux(int, const Point &x) const override {
return VectorField(-1.0, 0.0);
}
string Name() const override { return string("Stochastic Laplace 1D"); }
};
class StochasticLaplace2D : public StochasticProblem {
public:
explicit StochasticLaplace2D(StochasticFields &stochFields) :
StochasticProblem(stochFields) {} StochasticProblem(stochFields) {}
Scalar Solution(const Point &x) const override { return -x[1]; } Scalar Solution(const Point &x) const override { return -x[1]; }
...@@ -16,9 +30,11 @@ public: ...@@ -16,9 +30,11 @@ public:
StochasticProblem *getStochasticProblem(StochasticFields &stochFields, StochasticProblem *getStochasticProblem(StochasticFields &stochFields,
string StochProblem) { string StochProblem) {
if (StochProblem.empty()) if (StochProblem.empty()) ReadConfig(Settings, "Problem", StochProblem);
ReadConfig(Settings, "StochProblem", StochProblem); if (StochProblem == "StochasticLaplace1D")
if (StochProblem == "StochLaplace2D") return new StochasticProblem2D(stochFields); return new StochasticLaplace1D(stochFields);
if (StochProblem == "StochasticLaplace2D")
return new StochasticLaplace2D(stochFields);
else Exit("This StochProblem is not implemented") else Exit("This StochProblem is not implemented")
return nullptr; return nullptr;
} }
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