Commit 97e5e767 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactored DGTransportAssemble,

Equipped IStochasticTransportAssemble,
Implemented Missing functions in PDESolver (plot still missing),
Implemented GaussHat and StochasticGaussHat (not defined yet)
wrote tet
parent bf787c95
......@@ -282,3 +282,113 @@ void DGTransportAssemble::AssembleTransfer(TransferMatrix &TM) const {
TM(c(), c.Child(k))[0] = 1;
}
double DGTransportAssemble::Mass(const Vector &u) const {
Scalar e = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
e += w * U;
}
}
return PPM->Sum(e);
}
void DGTransportAssemble::SetInitialValue(Vector &u) const {
u = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
row r = u.find_row(c());
DGElement Elem(*disc, u, c);
double lambda = 1e-9;
for (int j = 0; j < Elem.NodalPoints(); ++j)
u(r, j) =
problem->Solution(0, lambda * c() + (1 - lambda) * Elem.NodalPoint(j));
}
Accumulate(u);
}
void DGTransportAssemble::PrintMatrixInfo(Matrix &A, int diagonal) const {
for (cell c = A.cells(); c != A.cells_end(); ++c) {
DGElement elem(*disc, A, c);
DGRowEntries A_c(A, c, c);
for (int i = 0; i < elem.NodalPoints(); ++i) {
for (int j = 0; j < elem.NodalPoints(); ++j)
mout << A_c(i, j) << " ";
mout << endl;
}
mout << endl;
if (diagonal) continue;
for (int f = 0; f < c.Faces(); ++f) {
cell cf = c;
if (A.GetMesh().onBndDG(c, f))
continue;
else
cf = A.GetMesh().find_neighbour_cell(c, f);
DGRowEntries A_cf(A, c, cf);
for (int i = 0; i < elem.NodalPoints(); ++i) {
for (int j = 0; j < elem.NodalPoints(); ++j)
mout << A_cf(i, j) << " ";
mout << endl;
}
mout << endl;
}
mout << "------------------------" << endl;
}
}
void DGTransportAssemble::VtkPlotting_cell(double t,
const Vector &u,
char *filename) const {
Vector u_tmp(u);
for (cell c = u_tmp.cells(); c != u_tmp.cells_end(); ++c) {
DGElement elem(*disc, u_tmp, c);
Scalar U = 0.0;
double a = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
U += w * elem.Value(q, u);
a += w;
}
U *= (1 / a);
row r = u.find_row(c());
u_tmp(r)[0] = U;
}
plot->celldata(u_tmp, 1);
NumberName("U", filename, step);
plot->vtk_celldata(filename, 0);
plot->gnu_celldata(filename, 0);
}
void DGTransportAssemble::VtkPlotting(double t, const Vector &u) const {
// Todo, refactor plot class should be able to handle everything
const Mesh &Mp = plot->GetMesh();
const Mesh &Mu = u.GetMesh();
if (Mp.CellCount() != Mu.CellCount()) return;
vout (2) << " => VtkPlotting result of step " << step << ".\n";
char filename[128];
VtkPlotting_cell(t, u, filename);
}
void DGTransportAssemble::PrintInfo() const {
mout.PrintInfo("Assemble", verbose,
PrintInfoEntry("Name", Name()),
PrintInfoEntry("Problem", problem->Name()),
PrintInfoEntry("Discretization", disc->DiscName()));
}
void DGTransportAssemble::FinishTimeStep(double t, Vector &u, bool SpecialTimeStep) {
step++;
t = t;
PrintInfo(u);
// VtkPlotting(t, u);
}
void DGTransportAssemble::Initialize(double t, Vector &u) {
u = 0;
SetInitialValue(u);
PrintInfo(u);
// VtkPlotting(t, u);
}
......@@ -31,23 +31,9 @@ public:
return disc;
}
void Initialize(double t, Vector &u) override {
u = 0;
SetInitialValue(u);
PrintInfo(u);
// VtkPlotting(t, u);
}
void SetGamma(Scalar gamma) {
this->gamma = gamma;
}
void Initialize(double t, Vector &u) override;
void FinishTimeStep(double t, Vector &u, bool SpecialTimeStep) override {
step++;
t = t;
PrintInfo(u);
// VtkPlotting(t, u);
}
void FinishTimeStep(double t, Vector &u, bool SpecialTimeStep) override;
const char *Name() const override { return "DGTransportAssemble"; }
......@@ -64,12 +50,7 @@ public:
mout << endl;
}
void PrintInfo() const {
mout.PrintInfo("Assemble", verbose,
PrintInfoEntry("Name", Name()),
PrintInfoEntry("Problem", problem->Name()),
PrintInfoEntry("Discretization", disc->DiscName()));
}
void PrintInfo() const override;
double Residual(const Vector &u, Vector &b) const override;
......@@ -85,95 +66,19 @@ public:
void Energy(const cell &c, const Vector &u, double &energy) const override;
double Mass(const Vector &u) const {
Scalar e = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
e += w * U;
}
}
return PPM->Sum(e);
}
double Mass(const Vector &u) const override;
double Error(double t, const Vector &u) const;
double Error(double t, const Vector &u) const override;
std::pair<double, double> InFlowOutFlowRate(const Vector &u) const;
RatePair InFlowOutFlowRate(const Vector &u) const override;
void SetInitialValue(Vector &u) const {
u = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
row r = u.find_row(c());
DGElement Elem(*disc, u, c);
double lambda = 1e-9;
for (int j = 0; j < Elem.NodalPoints(); ++j)
u(r, j) =
problem->Solution(0, lambda * c() + (1 - lambda) * Elem.NodalPoint(j));
}
Accumulate(u);
}
void SetInitialValue(Vector &u) const;
void VtkPlotting(double t, const Vector &u) const {
const Mesh &Mp = plot->GetMesh();
const Mesh &Mu = u.GetMesh();
if (Mp.CellCount() != Mu.CellCount()) return;
vout (2) << " => VtkPlotting result of step " << step << ".\n";
char filename[128];
VtkPlotting_cell(t, u, filename);
}
void VtkPlotting(double t, const Vector &u) const;
void VtkPlotting_cell(double t, const Vector &u, char *filename) const {
Vector u_tmp(u);
for (cell c = u_tmp.cells(); c != u_tmp.cells_end(); ++c) {
DGElement elem(*disc, u_tmp, c);
Scalar U = 0.0;
double a = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
U += w * elem.Value(q, u);
a += w;
}
U *= (1 / a);
row r = u.find_row(c());
u_tmp(r)[0] = U;
}
plot->celldata(u_tmp, 1);
NumberName("U", filename, step);
plot->vtk_celldata(filename, 0);
plot->gnu_celldata(filename, 0);
}
void VtkPlotting_cell(double t, const Vector &u, char *filename) const;
void PrintMatrixInfo(Matrix &A, int diagonal) const {
for (cell c = A.cells(); c != A.cells_end(); ++c) {
DGElement elem(*disc, A, c);
DGRowEntries A_c(A, c, c);
for (int i = 0; i < elem.NodalPoints(); ++i) {
for (int j = 0; j < elem.NodalPoints(); ++j)
mout << A_c(i, j) << " ";
mout << endl;
}
mout << endl;
if (diagonal) continue;
for (int f = 0; f < c.Faces(); ++f) {
cell cf = c;
if (A.GetMesh().onBndDG(c, f))
continue;
else
cf = A.GetMesh().find_neighbour_cell(c, f);
DGRowEntries A_cf(A, c, cf);
for (int i = 0; i < elem.NodalPoints(); ++i) {
for (int j = 0; j < elem.NodalPoints(); ++j)
mout << A_cf(i, j) << " ";
mout << endl;
}
mout << endl;
}
mout << "------------------------" << endl;
}
}
void PrintMatrixInfo(Matrix &A, int diagonal) const;
void SetExactSolution(double t, Vector &u_ex) const;
......
......@@ -6,13 +6,18 @@
#include "problem/StochasticTransportProblem.hpp"
typedef std::pair<double, double> RatePair;
class IStochasticTransportAssemble : public TAssemble {
protected:
InfoEntries entries;
IStochasticTransportProblem *problem;
public:
explicit IStochasticTransportAssemble(IStochasticTransportProblem *problem) :
problem(problem) {
InfoEntries problemEntries = problem->GetInfoEntries();
};
void DrawSample(const SampleID &id) {
......@@ -26,6 +31,14 @@ public:
virtual IDiscretization *GetDisc() = 0;
virtual ~IStochasticTransportAssemble() = default;
virtual double Energy(const Vector &u) const = 0;
virtual double Mass(const Vector &u) const = 0;
virtual double Error(double t, const Vector &u) const = 0;
virtual RatePair InFlowOutFlowRate(const Vector &u) const = 0;
};
#endif //ISTOCHASTICTRANSPORTASSEMBLE_HPP
......@@ -35,6 +35,14 @@ void TransportPDESolver::run(SampleSolution &solution) {
}
void TransportPDESolver::computeQ(SampleSolution &solution) {
if (quantity == "Energy") solution.Q = assemble->Energy(solution.U);
else if (quantity == "Mass") solution.Q = assemble->Mass(solution.U);
else if (quantity == "Error") solution.Q = assemble->Error(0, solution.U); // Todo
else if (quantity == "Inflow") solution.Q = assemble->InFlowOutFlowRate(solution.U).first;
else if (quantity == "Outflow") solution.Q = assemble->InFlowOutFlowRate(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) {
......
......@@ -178,8 +178,9 @@ public:
TransportPDESolver(IStochasticTransportAssemble *assemble) :
PDESolver(),
assemble(assemble),
// T(assemble) Todo make T and t0 Problem Properties
solver(Solver(GetPC("PointBlockGaussSeidel"), "GMRES")),
timeInt(TimeIntegrator(solver, 4)) {
timeInt(TimeIntegrator(solver, -102)) {
config.get("CFL", CFL);
config.get("t0", t0);
......
......@@ -15,7 +15,11 @@ CreateStochasticTransportProblem(std::string problemName, Meshes &meshes) {
return new StochasticPollutionCosHat2D(meshes);
if (problemName == "StochasticPollutionMollifiedBar2D")
return new StochasticPollutionMollifiedBar2D(meshes);
if (problemName == "StochasticInitialConditions2D")
return nullptr;
if (problemName == "GaussHat2D")
return new GaussHat2D(meshes);
if (problemName == "StochasticGaussHat2D")
return new StochasticGaussHat2D(meshes);
Exit(problemName + " not found")
}
......@@ -255,7 +255,69 @@ public:
string Name() const override { return "RiemannTransport1D"; }
};
class GaussHat2D : public IStochasticTransportProblem {
private:
double mu = 0.0;
double sigma = 1.0; // (sigma squared)
VectorField vectorField;
public:
GaussHat2D(Meshes &meshes) :
vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
IStochasticTransportProblem(GeneratorNames{"DummyGenerator"},
meshes) {}
Scalar Solution(double t, const Point &x) const override {
// Todo 2D Gaussfunction for t = 0
// 1 / (sigma * sqrt(2 * PI)) exp ()
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
return 1.0;
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return vectorField;
}
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return vectorField * N;
}
string Name() const override { return "GaussHat2D"; }
};
class StochasticGaussHat2D : public IStochasticTransportProblem {
private:
double mu = 0.0;
double sigma = 1.0; // (sigma squared)
VectorField vectorField;
public:
StochasticGaussHat2D(Meshes &meshes) :
vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
IStochasticTransportProblem(GeneratorNames{"DummyGenerator"},
meshes) {}
Scalar Solution(double t, const Point &x) const override {
// Todo 2D Gaussfunction for t = 0
// 1 / (sigma * sqrt(2 * PI)) exp ()
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
return 1.0;
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return vectorField;
}
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return vectorField * N;
}
string Name() const override { return "GaussHat2D"; }
};
//class StochasticInflowTransport1D : public IStochasticTransportProblem {
//public:
......
......@@ -2,7 +2,7 @@
#include "TestEnvironment.hpp"
constexpr double tolerance_pde_solver_test = 1e-10;
constexpr double PDESOLVER_TEST_TOLERANCE = 1e-10;
struct TestParams {
std::string model;
......@@ -36,6 +36,7 @@ protected:
dummyID.level.pLevel,
dummyID.level.fine);
pdeSolver = PDESolverCreator(*meshes).
WithModel(GetParam().model).
WithProblem(GetParam().problem).Create();
......@@ -50,9 +51,9 @@ protected:
}
};
class TestLaplace1D : public TestPDESolver {};
INSTANTIATE_TEST_CASE_P(TestPDESolver, TestPDESolver, Values(
// 1D Laplace
INSTANTIATE_TEST_CASE_P(TestPDESolver, TestLaplace1D, Values(
TestParams{"LagrangeElliptic", "Laplace1D", "L2Error", 0.0},
// TestParams{"LagrangeElliptic", "Laplace1D", "EnergyError", 0.0},
TestParams{"LagrangeElliptic", "Laplace1D", "FluxError", 0.0},
......@@ -63,15 +64,24 @@ INSTANTIATE_TEST_CASE_P(TestPDESolver, TestPDESolver, Values(
// TestParams{"MixedElliptic", "Laplace1D", "EnergyError", 0.0},
TestParams{"MixedElliptic", "Laplace1D", "FluxError", 0.0},
TestParams{"MixedElliptic", "Laplace1D", "Inflow", -1.0},
TestParams{"MixedElliptic", "Laplace1D", "Outflow", 1.0},
TestParams{"MixedElliptic", "Laplace1D", "Outflow", 1.0}
// TestParams{"HybridElliptic", "Laplace1D", "L2CellAverageError", 0.0},
// TestParams{"HybridElliptic", "Laplace1D", "Inflow", -1.0},
// TestParams{"HybridElliptic", "Laplace1D", "Outflow", 1.0},
// TestParams{"HybridElliptic", "Laplace1D", "EnergyError", 0.0},
// TestParams{"HybridElliptic", "Laplace1D", "FluxError", 0.0},
));
TEST_P(TestLaplace1D, TestRun) {
SampleSolution solution(*mGraphs, dummyID);
pdeSolver->Run(solution);
ASSERT_NEAR(solution.Q, GetParam().Q, PDESOLVER_TEST_TOLERANCE);
}
class TestLaplace2D : public TestPDESolver {};
// 2D Laplace
INSTANTIATE_TEST_CASE_P(TestPDESolver, TestLaplace2D, Values(
TestParams{"LagrangeElliptic", "Laplace2D", "L2Error", 0.0},
TestParams{"LagrangeElliptic", "Laplace2D", "EnergyError", 0.0},
TestParams{"LagrangeElliptic", "Laplace2D", "FluxError", 0.0},
......@@ -95,14 +105,25 @@ INSTANTIATE_TEST_CASE_P(TestPDESolver, TestPDESolver, Values(
TestParams{"DGElliptic", "Laplace2D", "FluxError", 0.0},
TestParams{"DGElliptic", "Laplace2D", "Inflow", -1.0},
TestParams{"DGElliptic", "Laplace2D", "Outflow", 1.0}
));
TEST_P(TestLaplace2D, TestRun) {
SampleSolution solution(*mGraphs, dummyID);
pdeSolver->Run(solution);
ASSERT_NEAR(solution.Q, GetParam().Q, PDESOLVER_TEST_TOLERANCE);
}
class TestGaussHat : public TestPDESolver {};
// TestParams{"DGTransport", "CircleWave2D", "Mass", 0.0}
INSTANTIATE_TEST_CASE_P(TestPDESolver, TestGaussHat, Values(
TestParams{"DGTransport", "GaussHat2D", "Mass", 1.0}
));
TEST_P(TestPDESolver, TestRun) {
TEST_P(TestGaussHat, TestRun) {
SampleSolution solution(*mGraphs, dummyID);
pdeSolver->Run(solution);
ASSERT_NEAR(solution.Q, GetParam().Q, tolerance_pde_solver_test);
ASSERT_NEAR(solution.Q, GetParam().Q, PDESOLVER_TEST_TOLERANCE);
// Todo ASSERT_NEAR(solution.Cost, GetParam().Cost, PDESOLVER_TEST_TOLERANCE);
}
int main(int argc, char **argv) {
......
Supports Markdown
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