Commit e4cd5835 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Merge branch '51-implement-and-test-parabolic-pde-solver' into 'feature'

Resolve "Implement and test Parabolic PDE Solver"

Closes #51

See merge request !96
parents bc5f0a40 6e4252ad
Pipeline #179952 passed with stages
in 42 minutes
......@@ -2,8 +2,7 @@
#include "Plotting.hpp"
void DGReactionAssemble::Dirichlet(double t, Vector &u) const {
// Todo has to be checked
void DGReactionAssemble::Initialize(Vector &u) const {
u.ClearDirichletFlags();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
RowBndValues u_c(u, c);
......@@ -17,7 +16,7 @@ void DGReactionAssemble::Dirichlet(double t, Vector &u) const {
if (F > -1e-6) continue;
for (int j = 0; j < u.NodalPointsOnFace(c, i); ++j) {
int k = u.NodalPointOnFace(c, i, j);
u_c(k) = problem->Concentration(t, E.NodalPoint(k));
u_c(k) = problem->Concentration(Time(), E.NodalPoint(k));
u_c.D(k) = true;
}
}
......@@ -25,9 +24,9 @@ void DGReactionAssemble::Dirichlet(double t, Vector &u) const {
u.DirichletConsistent();
}
double DGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const {
double DGReactionAssemble::Residual(const Vector &u, Vector &r) const {
r = 0;
const Vector &u_old = U_old();
const Vector &u_old = PreviousSolution();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
int sd = c.Subdomain();
RowBndValues r_c(r, c);
......@@ -48,10 +47,8 @@ double DGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const
for (int i = 0; i < elem.NodalPoints(); ++i) {
Scalar U_i = elem.Value(q, i);
VectorField DU_i = elem.Derivative(q, i);
r(c(), i) += w * (K_DU * DU_i
+ (B * DU) * U_i
+ (1.0 / dt_) * (U - Uold) * U_i
- (R + f) * U_i);
r(c(), i) += w * (K_DU * DU_i + (B * DU) * U_i
+ (1.0 / StepSize()) * (U - Uold) * U_i - (R + f) * U_i);
}
}
for (int f = 0; f < c.Faces(); ++f) {
......@@ -64,7 +61,7 @@ double DGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const
Point N = felem.QNormal(q);
Scalar BN = problem->FaceConvection(c, f, N, z);
Tensor K = problem->Diffusion(z);
Scalar U_diff = felem.Value(q, u) - problem->Concentration(t, z);
Scalar U_diff = felem.Value(q, u) - problem->Concentration(Time(), z);
VectorField K_DU = K * felem.Derivative(q, u);
if (BN < 0) {
Scalar F = K_DU * N;
......@@ -155,7 +152,7 @@ void DGReactionAssemble::Jacobi(const Vector &u, Matrix &A) const {
VectorField DU_j = elem.Derivative(q, j);
A_c(i, j) += w * (K_DU_i * DU_j
+ (B * DU_j) * U_i
+ (1.0 / dt_) * U_i * U_j
+ (1.0 / StepSize()) * U_i * U_j
- DR * U_i * U_j);
}
}
......
......@@ -44,9 +44,9 @@ public:
PrintInfoEntry("Discretization", disc.DiscName()));
}
void Dirichlet(double t, Vector &u) const override;
void Initialize(Vector &u) const override;
double Residual(double t, const Vector &u, Vector &r) const override;
double Residual(const Vector &u, Vector &r) const override;
void Jacobi(const Vector &u, Matrix &A) const override;
......
......@@ -36,51 +36,31 @@ public:
problem->GetStepSize(u.GetMesh().MaxMeshWidth()));
}
virtual void Dirichlet(double t, Vector &u) const override = 0;
virtual double Energy(const Vector &u) const override { return 0.0; };
virtual double Residual(double t, const Vector &u, Vector &r) const override = 0;
double Energy(const Vector &u) const override { return 0.0; };
virtual double Mass(const Vector &u) const = 0;
virtual RatePair InflowOutflow(const Vector &u) const = 0;
void PrintIteration(const Vector &u) const override {
if (verbose == 0 || Step() % verbose != 0) return;
std::pair<double, double> rate = InflowOutflow(u);
mout << "n=" << std::setw(3) << std::left << step
<< " t=" << std::setw(13) << std::left << t_
<< " 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;
// if (verbose == 0 || Step() % verbose != 0) return;
// std::pair<double, double> rate = InflowOutflow(u);
// mout.PrintIteration(verbose,
// PrintIterEntry("n", Step(), int(log10(MaxStep())), 1),
// PrintIterEntry("t", Time(), 13, 1),
// PrintIterEntry("Mass", Mass(u), 13, 1),
// PrintIterEntry("IFR", rate.first, 13, 1),
// PrintIterEntry("OFR", rate.second, 13, 1)
// );
mout.PrintIteration(
verbose,
PrintIterEntry("n", Step(), int(log10(MaxStep())), 1),
PrintIterEntry("t", Time(), 13, 1),
PrintIterEntry("Mass", Mass(u), 13, 1),
PrintIterEntry("IFR", rate.first, 13, 1),
PrintIterEntry("OFR", rate.second, 13, 1)
);
};
void Initialize(double t, Vector &u) override {
SetInitialValue(u);
}
void FinishTimeStep(double t, Vector &u) override {
++step;
PlotIteration(u);
PrintIteration(u);
}
const char *Name() const override = 0;
virtual void PrintInfo() const = 0;
};
#endif //ISTOCHASTICREACTIONASSEMBLE_HPP
......@@ -2,7 +2,7 @@
#include "Plotting.hpp"
void PGReactionAssemble::Dirichlet(double t, Vector &u) const {
void PGReactionAssemble::Initialize(Vector &u) const {
u.ClearDirichletFlags();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
RowBndValues u_c(u, c);
......@@ -12,12 +12,11 @@ void PGReactionAssemble::Dirichlet(double t, Vector &u) const {
for (int i = 0; i < c.Faces(); ++i) {
if (u_c.bc(i) == -1) continue;
ScalarFaceElement FE(u, c, i);
Scalar F = problem->FaceConvection(c, i, FE.QNormal(0),
FE.QPoint(0));
Scalar F = problem->FaceConvection(c, i, FE.QNormal(0), FE.QPoint(0));
if (F > -1e-6) continue;
for (int j = 0; j < u.NodalPointsOnFace(c, i); ++j) {
int k = u.NodalPointOnFace(c, i, j);
u_c(k) = problem->Concentration(t, E.NodalPoint(k));
u_c(k) = problem->Concentration(Time(), E.NodalPoint(k));
u_c.D(k) = true;
}
}
......@@ -25,9 +24,8 @@ void PGReactionAssemble::Dirichlet(double t, Vector &u) const {
u.DirichletConsistent();
}
double PGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const {
double PGReactionAssemble::Residual(const Vector &u, Vector &r) const {
r = 0;
const Vector &u_old = U_old();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
ScalarElement E(u, c);
RowBndValues r_c(r, c);
......@@ -35,34 +33,28 @@ double PGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
Scalar U = E.Value(q, u);
Scalar Uold = E.Value(q, u_old);
Scalar Uold = E.Value(q, PreviousSolution());
Tensor K = problem->Diffusion(E.QPoint(q));
Scalar f = problem->Load(E.QPoint(q));
VectorField C = problem->CellConvection(
c, E.QPoint(q));
VectorField C = problem->CellConvection(c, E.QPoint(q));
Scalar R = problem->Reaction(E.QPoint(q), U);
VectorField DU = E.Derivative(q, u);
for (unsigned int i = 0; i < E.size(); ++i) {
Scalar U_i = E.Value(q, i);
VectorField DU_i = E.Derivative(q, i);
r_c(i) += w * ((K * DU) * DU_i + U_i * (C * DU)
+ (1.0 / dt_) * (U - Uold) * U_i
- (R + f) * U_i
+ delta_h *
((C * DU) + (1.0 / dt_) * (U - Uold) - (R + f))
* (C * DU_i));
+ (1.0 / StepSize()) * (U - Uold) * U_i - (R + f) * U_i
+ delta_h * ((C * DU) + (1.0 / StepSize()) * (U - Uold) - (R + f)) * (C * DU_i));
}
BFParts bnd(u.GetMesh(), c);
if (bnd.onBnd()) continue;
int sd = c.Subdomain();
for (int f = 0; f < c.Faces(); ++f) {
if (bnd[f] == -1) continue;
ScalarFaceElement FE(u, c, f);
RowValues r_f(r, FE);
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Scalar F = problem->FaceConvection(c, f, FE.QNormal(q),
FE.QPoint(q));
Scalar F = problem->FaceConvection(c, f, FE.QNormal(q), FE.QPoint(q));
if (F < 1e-6) continue;
Scalar U = 0;
for (unsigned int i = 0; i < FE.size(); ++i) {
......@@ -101,8 +93,8 @@ void PGReactionAssemble::Jacobi(const Vector &u, Matrix &A) const {
VectorField DU_j = E.Derivative(q, j);
Scalar U_j = E.Value(q, j);
A_c(i, j) += w * ((K * DU_i) * DU_j + (C * DU_j) * U_i
+ (1.0 / dt_) * U_i * U_j - DR * U_i * U_j
+ delta_h * ((C * DU_j) + (1.0 / dt_) * U_j - DR * U_j) * (C * DU_i));
+ (1.0 / StepSize()) * U_i * U_j - DR * U_i * U_j
+ delta_h * ((C * DU_j) + (1.0 / StepSize()) * U_j - DR * U_j) * (C * DU_i));
}
}
BFParts bnd(u.GetMesh(), c);
......
......@@ -36,9 +36,9 @@ public:
PrintInfoEntry("Discretization", disc.DiscName()));
}
void Dirichlet(double t, Vector &u) const override;
void Initialize(Vector &u) const override;
double Residual(double t, const Vector &u, Vector &r) const override;
double Residual(const Vector &u, Vector &r) const override;
void Jacobi(const Vector &u, Matrix &A) const override;
......
Subproject commit f3907f7be890b4599c2fe9fab7f4bf8b8f421d63
Subproject commit c61e4cd7786f5279f914ae730e4232f80eb06915
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