Skip to content
GitLab
Commit c8aab1f4 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

adaption to 209

parent ccbd5a27
Hide whitespace changes
Inline Side-by-side
mlmc/src/pdesolver/IStochasticTransportAssemble.hpp
View file @ c8aab1f4
......@@ -11,7 +11,7 @@
typedef std::pair<double, double> RatePair;
class IStochasticTransportAssemble : public TAssemble {
class IStochasticLinearTransportAssemble : public ILinearTimeAssemble {
protected:
Plot *plot;
......@@ -20,7 +20,8 @@ protected:
IStochasticTransportProblem *problem;
public:
explicit IStochasticTransportAssemble(IStochasticTransportProblem *problem) :
explicit IStochasticLinearTransportAssemble(IStochasticTransportProblem *problem) :
ILinearTimeAssemble(TimeSeries()),
problem(problem) {
};
......@@ -44,7 +45,7 @@ public:
virtual IDiscretization *GetDisc() = 0;
virtual ~IStochasticTransportAssemble() = default;
virtual ~IStochasticLinearTransportAssemble() = default;
virtual double Energy(const Vector &u) const = 0;
......@@ -55,4 +56,9 @@ public:
virtual RatePair InFlowOutFlowRate(const Vector &u) const = 0;
};
class IStochasticNonLinearTransportAssemble : public INonLinearTimeAssemble {
};
#endif //ISTOCHASTICTRANSPORTASSEMBLE_HPP
mlmc/src/pdesolver/PDESolver.cpp
View file @ c8aab1f4
......@@ -98,7 +98,7 @@ void EllipticPDESolver::createOtherMatrixGraphs() {
void TransportPDESolver::run(SampleSolution &solution) {
auto timeSeries = assemble->GetTimeSeries(solution.U);
timeInt.operator()(assemble, timeSeries, solution.U);
timeInt->Method(assemble, solution.U);
}
void TransportPDESolver::computeQ(SampleSolution &solution) {
......
mlmc/src/pdesolver/PDESolver.hpp
View file @ c8aab1f4
......@@ -206,11 +206,9 @@ public:
class TransportPDESolver : public PDESolver {
private:
Solver solver;
TimeIntegrator timeInt;
TimeIntegrator *timeInt;
IStochasticTransportAssemble *assemble;
IStochasticLinearTransportAssemble *assemble;
protected:
void run(SampleSolution &solution) override;
......@@ -224,25 +222,28 @@ protected:
void setUpPlot(SampleSolution &solution) override;
public:
TransportPDESolver(IStochasticTransportAssemble *assemble,
TransportPDESolver(IStochasticLinearTransportAssemble *assemble,
Meshes &meshes,
const std::string &quantity = "L2",
const std::string &costMeasure = "size") :
PDESolver(meshes, quantity, costMeasure), assemble(assemble),
solver(Solver(GetPC("PointBlockJacobi_dG"), "GMRES")),
timeInt(TimeIntegrator(solver, -102)) {
timeInt(TimeIntegratorCreator(IMPLICIT_MIDPOINT).
WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))).
WithLinearSolver(new GMRES(GetPC("PointBlockJacobi_dG"))).
Create()) {
solMGraphs = CreateSolutionMatrixGraphs();
solver.PrintInfo();
PrintInfo();
}
~TransportPDESolver() {
if (!assemble) delete assemble;
if (!solMGraphs) delete solMGraphs;
delete timeInt;
}
IAssemble *GetAssemble() const override {
return assemble;
return nullptr;
// return assemble;
}
IDiscretization *GetDisc() const override {
......
mlmc/src/pdesolver/PDESolverCreator.cpp
View file @ c8aab1f4
......@@ -36,7 +36,7 @@ PDESolver *PDESolverCreator::Create(Meshes &meshes) {
if (_model == "DGTransport")
return new TransportPDESolver(
new DGTransportAssemble(
new DGLinearTransportAssemble(
new DGDiscretization(meshes, _degree),
CreateStochasticTransportProblem(_problem, meshes)
), meshes, _quantity, _costMeasure
......
mlmc/src/pdesolver/assembling/DGTransportAssemble.cpp
View file @ c8aab1f4
#include "DGTransportAssemble.hpp"
double DGTransportAssemble::Residual(const Vector &u, Vector &b) const {
b = 0;
auto *massMatrix = new Matrix(u);
MassMatrix(*massMatrix);
auto *fluxMatrix = new Matrix(u);
SystemMatrix(*fluxMatrix);
Vector fluxMatrixU(u);
fluxMatrixU = *fluxMatrix * u;
fluxMatrixU *= -dt_;
b = (*massMatrix * u + fluxMatrixU);
b -= *massMatrix * U_old();
Vector rhs(b);
RHS(t_, rhs);
b -= dt_ * rhs;
b.ClearDirichletValues();
Collect(b);
delete fluxMatrix;
delete massMatrix;
return b.norm();
}
void DGTransportAssemble::Jacobi(const Vector &u, Matrix &A) const {
Matrix *massMatrix = new Matrix(u);
MassMatrix(*massMatrix);
Matrix *fluxMatrix = new Matrix(u);
SystemMatrix(*fluxMatrix);
A = *massMatrix;
A += -dt_ * (*fluxMatrix);
delete fluxMatrix;
delete massMatrix;
}
void DGTransportAssemble::MassMatrix(Matrix &massMatrix) const {
void DGLinearTransportAssemble::MassMatrix(Matrix &massMatrix) const {
massMatrix = 0;
for (cell c = massMatrix.cells(); c != massMatrix.cells_end(); ++c) {
DGElement elem(*disc, massMatrix, c);
......@@ -57,7 +19,7 @@ void DGTransportAssemble::MassMatrix(Matrix &massMatrix) const {
}
}
void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
void DGLinearTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
systemMatrix = 0;
for (cell c = systemMatrix.cells(); c != systemMatrix.cells_end(); ++c) {
DGElement elem(*disc, systemMatrix, c);
......@@ -167,7 +129,7 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
}
}
void DGTransportAssemble::RHS(double t, Vector &rhs) const {
void DGLinearTransportAssemble::RHS(double t, Vector &rhs) const {
rhs = 0;
if (!problem->RHS()) return;
for (cell c = rhs.cells(); c != rhs.cells_end(); ++c) {
......@@ -191,7 +153,7 @@ void DGTransportAssemble::RHS(double t, Vector &rhs) const {
}
}
double DGTransportAssemble::Energy(const Vector &u) const {
double DGLinearTransportAssemble::Energy(const Vector &u) const {
double energy = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
......@@ -204,7 +166,7 @@ double DGTransportAssemble::Energy(const Vector &u) const {
return 0.5 * PPM->Sum(energy);
}
void DGTransportAssemble::Energy(const cell &c, const Vector &u, double &energy) const {
void DGNonLinearTransportAssemble::Energy(const cell &c, const Vector &u, double &energy) const {
DGElement elem(*disc, u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
......@@ -213,7 +175,7 @@ void DGTransportAssemble::Energy(const cell &c, const Vector &u, double &energy)
}
}
double DGTransportAssemble::Error(double t, const Vector &u) const {
double DGLinearTransportAssemble::Error(double t, const Vector &u) const {
double err = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
......@@ -227,7 +189,7 @@ double DGTransportAssemble::Error(double t, const Vector &u) const {
return sqrt(PPM->Sum(err));
}
std::pair<double, double> DGTransportAssemble::InFlowOutFlowRate(const Vector &u) const {
std::pair<double, double> DGLinearTransportAssemble::InFlowOutFlowRate(const Vector &u) const {
double inflow = 0.0;
double outflow = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
......@@ -263,7 +225,7 @@ std::pair<double, double> DGTransportAssemble::InFlowOutFlowRate(const Vector &u
return {inflow, outflow};
}
void DGTransportAssemble::SetExactSolution(double t, Vector &u_ex) const {
void DGLinearTransportAssemble::SetExactSolution(double t, Vector &u_ex) const {
u_ex = 0;
for (cell c = u_ex.cells(); c != u_ex.cells_end(); ++c) {
row r = u_ex.find_row(c());
......@@ -274,7 +236,7 @@ void DGTransportAssemble::SetExactSolution(double t, Vector &u_ex) const {
Accumulate(u_ex);
}
void DGTransportAssemble::AssembleTransfer(TransferMatrix &TM) const {
void DGNonLinearTransportAssemble::AssembleTransfer(TransferMatrix &TM) const {
TM = 0;
const matrixgraph &cg = TM.CoarseMatrixGraph();
for (cell c = cg.cells(); c != cg.cells_end(); ++c)
......@@ -282,7 +244,7 @@ void DGTransportAssemble::AssembleTransfer(TransferMatrix &TM) const {
TM(c(), c.Child(k))[0] = 1;
}
double DGTransportAssemble::Mass(const Vector &u) const {
double DGLinearTransportAssemble::Mass(const Vector &u) const {
Scalar e = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
......@@ -295,7 +257,7 @@ double DGTransportAssemble::Mass(const Vector &u) const {
return PPM->Sum(e);
}
void DGTransportAssemble::SetInitialValue(Vector &u) const {
void DGLinearTransportAssemble::SetInitialValue(Vector &u) const {
u = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
row r = u.find_row(c());
......@@ -308,7 +270,7 @@ void DGTransportAssemble::SetInitialValue(Vector &u) const {
Accumulate(u);
}
void DGTransportAssemble::PrintMatrixInfo(Matrix &A, int diagonal) const {
void DGLinearTransportAssemble::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);
......@@ -338,9 +300,9 @@ void DGTransportAssemble::PrintMatrixInfo(Matrix &A, int diagonal) const {
}
}
void DGTransportAssemble::VtkPlotting_cell(double t,
const Vector &u,
char *filename) const {
void DGLinearTransportAssemble::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);
......@@ -361,7 +323,7 @@ void DGTransportAssemble::VtkPlotting_cell(double t,
plot->gnu_celldata(filename, 0);
}
void DGTransportAssemble::VtkPlotting(double t, const Vector &u) const {
void DGLinearTransportAssemble::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();
......@@ -371,24 +333,64 @@ void DGTransportAssemble::VtkPlotting(double t, const Vector &u) const {
VtkPlotting_cell(t, u, filename);
}
void DGTransportAssemble::PrintInfo() const {
void DGLinearTransportAssemble::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) {
void DGLinearTransportAssemble::FinishTimeStep(double t, Vector &u) {
step++;
t = t;
PrintInfo(u);
VtkPlotting(t, u);
}
void DGTransportAssemble::Initialize(double t, Vector &u) {
void DGLinearTransportAssemble::Initialize(Vector &u) {
u = 0;
SetInitialValue(u);
PrintInfo(u);
VtkPlotting(t, u);
VtkPlotting(timeSeries.FirstTStep(), u);
}
double DGNonLinearTransportAssemble::Residual(const Vector &u, Vector &b) const {
b = 0;
auto *massMatrix = new Matrix(u);
MassMatrix(*massMatrix);
auto *fluxMatrix = new Matrix(u);
SystemMatrix(*fluxMatrix);
Vector fluxMatrixU(u);
fluxMatrixU = *fluxMatrix * u;
fluxMatrixU *= -dt_;
b = (*massMatrix * u + fluxMatrixU);
b -= *massMatrix * U_old();
Vector rhs(b);
RHS(t_, rhs);
b -= dt_ * rhs;
b.ClearDirichletValues();
Collect(b);
delete fluxMatrix;
delete massMatrix;
return b.norm();
}
void DGNonLinearTransportAssemble::Jacobi(const Vector &u, Matrix &A) const {
Matrix *massMatrix = new Matrix(u);
MassMatrix(*massMatrix);
Matrix *fluxMatrix = new Matrix(u);
SystemMatrix(*fluxMatrix);
A = *massMatrix;
A += -dt_ * (*fluxMatrix);
delete fluxMatrix;
delete massMatrix;
}
mlmc/src/pdesolver/assembling/DGTransportAssemble.hpp
View file @ c8aab1f4
......@@ -9,15 +9,16 @@
#include <iomanip>
class DGTransportAssemble : public IStochasticTransportAssemble {
class DGLinearTransportAssemble : public IStochasticLinearTransportAssemble {
public:
DGDiscretization *disc;
double flux_alpha = 1.0;
double diffusion = 0.0;
DGTransportAssemble(IDiscretizationT<> *disc, IStochasticTransportProblem *problem)
: IStochasticTransportAssemble(problem),
DGLinearTransportAssemble(IDiscretizationT<> *disc,
IStochasticTransportProblem *problem)
: IStochasticLinearTransportAssemble(problem),
disc(dynamic_cast<DGDiscretization *> (disc)) {
config.get("flux_alpha", flux_alpha);
PrintInfo();
......@@ -27,10 +28,6 @@ public:
return disc;
}
void Initialize(double t, Vector &u) override;
void FinishTimeStep(double t, Vector &u, bool SpecialTimeStep) override;
const char *Name() const override { return "DGTransportAssemble"; }
void PrintInfo(const Vector &u) const override {
......@@ -46,11 +43,9 @@ public:
mout << endl;
}
void PrintInfo() const override;
double Residual(const Vector &u, Vector &b) const override;
void FinishTimeStep(double t, Vector &u) override;
void Jacobi(const Vector &u, Matrix &A) const override;
void PrintInfo() const override;
void MassMatrix(Matrix &massMatrix) const override;
......@@ -60,14 +55,14 @@ public:
double Energy(const Vector &u) const override;
void Energy(const cell &c, const Vector &u, double &energy) const override;
double Mass(const Vector &u) const override;
double Error(double t, const Vector &u) const override;
RatePair InFlowOutFlowRate(const Vector &u) const override;
void Initialize(Vector &u) override;
void SetInitialValue(Vector &u) const;
void VtkPlotting(double t, const Vector &u) const;
......@@ -77,115 +72,34 @@ public:
void PrintMatrixInfo(Matrix &A, int diagonal) const;
void SetExactSolution(double t, Vector &u_ex) const;
};
class DGNonLinearTransportAssemble : public IStochasticNonLinearTransportAssemble {
public:
IStochasticTransportProblem *problem;
DGDiscretization *disc;
Plot *plot;
double flux_alpha = 1.0;
double diffusion = 0.0;
DGNonLinearTransportAssemble(IDiscretizationT<> *disc,
IStochasticTransportProblem *problem)
: problem(problem), disc(dynamic_cast<DGDiscretization *> (disc)) {
config.get("flux_alpha", flux_alpha);
}
void Initialize(double t, Vector &u) override;
double Residual(const Vector &u, Vector &b) const override;
void Jacobi(const Vector &u, Matrix &A) const override;
double Energy(const Vector &u) const override;
void Energy(const cell &c, const Vector &u, double &energy) const override;
void AssembleTransfer(TransferMatrix &TM) const override;
};
#endif //TUTORIAL_DGTRANSPORT_HPP
//class DGTransportAssemble : public IStochasticTransportAssemble {
//protected:
// DGDiscretization *disc;
//
// double flux_alpha = 1.0;
// double diffusion = 0.0;
//
//public:
// DGTransportAssemble(IDiscretization *disc, IStochasticTransportProblem *problem) :
// IStochasticTransportAssemble(problem),
// disc(dynamic_cast<DGDiscretization *>(disc)) {
// config.get("flux_alpha", flux_alpha);
// config.get("diffusion", diffusion);
// }
//
// IDiscretization *GetDisc() override {
// return disc;
// };
//
// void PrintInfo() const {
// mout.PrintInfo("Assemble", 1,
// PrintInfoEntry("Name", Name()),
//// PrintInfoEntry("Problem", problem->Name()),
// PrintInfoEntry("Discretization", disc->DiscName()));
// }
//
// void PrintInfo(const Vector &u) const override {
//// std::pair<double, double> rate = InFlowOutFlowRate(u);
//// mout << "n=" << std::setw(3) << std::left << step
//// << " t=" << std::setw(7) << std::left << t_
//// << " Energy=" << std::setw(13) << std::left << Energy(u)
//// << " Mass=" << std::setw(13) << std::left << Mass(u);
//// if (abs(rate.first) >= 10e-10)
//// mout << " IFR=" << std::setw(13) << std::left << rate.first;
//// if (abs(rate.second) >= 10e-10)
//// mout << " OFR=" << std::setw(13) << std::left << rate.second;
//// mout << endl;
// }
//
// const char *Name() const override;
//
// double Residual(const Vector &u, Vector &b) const override;
//
// void Jacobi(const Vector &u, Matrix &A) const override;
//
// void MassMatrix(Matrix &massMatrix) const override;
//
// void SystemMatrix(Matrix &systemMatrix) const override;
//
// void RHS(double t, Vector &rhs) const override;
//
// double Energy(const Vector &u) const override;
//
// void Energy(const cell &c, const Vector &u, double &energy) const override;
//
// double Error(double t, const Vector &u) const;
//
// std::pair<double, double> InFlowOutFlowRate(const Vector &u) const;
//
// 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 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 SetExactSolution(double t, Vector &u_ex) const;
//
// void AssembleTransfer(TransferMatrix &TM) const override;
//};
//
//#endif //TUTORIAL_DGTRANSPORT_HPP
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

Imprint