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

updated reference on mpp submodule

parent 41a05672
Hide whitespace changes
Inline Side-by-side
mlmc/src/assemble/DGEllipticAssemble.C 0 → 100644
View file @ 347756ff
#include "DGEllipticAssemble.h"
DGEllipticAssemble::DGEllipticAssemble(const DGDiscretization &disc,
const ScalarProblem &problem)
: EllipticAssemble(disc, problem), disc(disc), problem(problem) {
ReadConfig(Settings, "penalty", penalty);
ReadConfig(Settings, "sign", sign);
}
DGEllipticAssemble::DGEllipticAssemble(const DGDiscretization &disc,
const ScalarProblem &problem, double penalty,
int sign)
: EllipticAssemble(disc, problem), disc(disc), problem(problem),
penalty(penalty), sign(sign) {}
Element *DGEllipticAssemble::getElement(const cell &c, const Vector &u) const {
return new DGElement(disc, u, c);
}
FaceElement *DGEllipticAssemble::getFaceElement(const cell &c, const Vector &u,
int face) const {
return new DGFaceElement(disc, u, c, face);
}
Scalar DGEllipticAssemble::getValue(const cell &c, const Vector &u,
Element *elem, int q) const {
auto dgEllipticElem = dynamic_cast<DGElement *>(elem);
return dgEllipticElem->Value(q, u);
}
Scalar DGEllipticAssemble::getFaceValue(const cell &c, const Vector &u,
FaceElement *faceElem, int q, int face) const {
auto dgEllipticFaceElem = dynamic_cast<DGFaceElement *>(faceElem);
return dgEllipticFaceElem->Value(q, u);
}
VectorField DGEllipticAssemble::getFlux(const cell &c, const Vector &u,
Element *elem, int q) const {
auto dgEllipticElem = dynamic_cast<DGElement *>(elem);
Tensor K = problem.Permeability(c.Subdomain(), elem->QPoint(q));
return K * dgEllipticElem->Derivative(q, u);
}
VectorField DGEllipticAssemble::getFaceFlux(const cell &c, const Vector &u,
FaceElement *faceElem, int q,
int face) const {
auto dgEllipticFaceElem = dynamic_cast<DGFaceElement *>(faceElem);
Tensor K = problem.Permeability(c.Subdomain(), dgEllipticFaceElem->QPoint(q));
return K * dgEllipticFaceElem->Derivative(q, u);
}
VectorField DGEllipticAssemble::getGradient(const cell &c, const Vector &u,
Element *elem, int q) const {
auto dgEllipticElem = dynamic_cast<DGElement *>(elem);
return dgEllipticElem->Derivative(q, u);
}
const char *DGEllipticAssemble::Name() const { return "Discontinuous-Galerkin"; }
void DGEllipticAssemble::Dirichlet(const cell &c, Vector &u) const {}
void DGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) const {
DGElement elem(disc, u, c);
RowBndValues r_c(r, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
const Point &z = elem.QPoint(q);
Tensor K = problem.Permeability(c.Subdomain(), z);
Scalar f = problem.Load(c.Subdomain(), z);
VectorField K_DU = K * elem.Derivative(q, u);
for (int i = 0; i < elem.dimension(); ++i) {
Scalar U_i = elem.Value(q, i);
VectorField DU_i = elem.Derivative(q, i);
r(c(), i) += w * (K_DU * DU_i - f * U_i);
}
}
for (int f = 0; f < c.Faces(); ++f) {
Scalar s = penalty / c.FaceArea(f);
DGFaceElement felem(disc, u, c, f);
if (checkIfBND(r, c, f)) {
int bnd = r_c.bc(f);
if (bnd == 2) {
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
const Point &z = felem.QPoint(q);
Scalar F = problem.Flux(bnd, z) * felem.QNormal(q);
for (int i = 0; i < felem.dimension(); ++i) {
Scalar U_i = felem.Value(q, i);
r(c(), i) -= w * U_i * F;
}
}
} else if (bnd == 1) {
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
const Point &z = felem.QPoint(q);
const Point &N = felem.QNormal(q);
Tensor K = problem.Permeability(c.Subdomain(), c());
Scalar U_diff = felem.Value(q, u) - problem.Solution(z);
VectorField K_DU = K * felem.Derivative(q, u);
Scalar F = K_DU * N;
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
r(c(), i) -= w * (F * phi_i +
U_diff * (NDphi_i * sign - s * phi_i));
}
}
}
} else {
cell cf = findNBCell(r, c, f);
if (cf() < c()) continue;
int f1 = findNBFaceID(r, c.Face(f), cf);
DGFaceElement felem_1(disc, r, cf, f1);
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
const Point &N = felem.QNormal(q);
Tensor K = problem.Permeability(c.Subdomain(), c());
Scalar U = felem.Value(q, u);
VectorField K_DU = K * felem.Derivative(q, u);
const Point &Qf_c = felem.QPoint(q);
int q1 = findQPointID(felem, felem_1, Qf_c);
Tensor K_1 = problem.Permeability(cf.Subdomain(), cf());
Scalar U_1 = felem_1.Value(q1, u);
VectorField K_DU_1 = K_1 * felem_1.Derivative(q1, u);
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
r(c(), i) += w * ((s * U - 0.5 * (K_DU * N)) * phi_i
- 0.5 * U * NDphi_i * sign);
r(c(), i) += w * ((-s * U_1 - 0.5 * (K_DU_1 * N)) * phi_i
+ 0.5 * U_1 * NDphi_i * sign);
}
for (int j = 0; j < felem_1.dimension(); ++j) {
Scalar phi_j = felem_1.Value(q1, j);
Scalar NDphi_j = (K_1 * felem_1.Derivative(q1, j)) * N;
r(cf(), j) += w * ((0.5 * K_DU_1 * N + s * U_1) * phi_j
+ 0.5 * U_1 * NDphi_j * sign);
r(cf(), j) += w * ((0.5 * K_DU * N - s * U) * phi_j
- 0.5 * U * NDphi_j * sign);
}
}
}
}
}
void DGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const {
DGElement elem(disc, A, c);
DGRowEntries A_c(A, c, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem.Permeability(c.Subdomain(), z);
for (int i = 0; i < elem.dimension(); ++i) {
VectorField K_DU_i = K * elem.Derivative(q, i);
for (int j = 0; j < elem.dimension(); ++j) {
VectorField DU_j = elem.Derivative(q, j);
A_c(i, j) += w * K_DU_i * DU_j;
}
}
}
BFParts bnd(u.GetMesh(), c);
for (int f = 0; f < c.Faces(); ++f) {
Scalar s = penalty / c.FaceArea(f);
DGFaceElement felem(disc, u, c, f);
if (checkIfBND(u, c, f)) {
if (bnd[f] == 1) {
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
const Point &z = felem.QPoint(q);
const Point &N = felem.QNormal(q);
Tensor K = problem.Permeability(c.Subdomain(), c());
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
for (int j = 0; j < felem.dimension(); ++j) {
Scalar phi_j = felem.Value(q, j);
Scalar NDphi_j = (K * felem.Derivative(q, j)) * N;
A_c(i, j) -= w * (NDphi_i * phi_j * sign +
phi_i * (NDphi_j - s * phi_j));
}
}
}
}
} else {
cell cf = findNBCell(u, c, f);
if (cf() < c()) continue;
DGRowEntries A_cf(A, c, cf);
DGRowEntries A_fc(A, cf, c);
DGRowEntries A_ff(A, cf, cf);
int f1 = findNBFaceID(u, c.Face(f), cf);
DGFaceElement felem_1(disc, u, cf, f1);
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
const Point &N = felem.QNormal(q);
Tensor K = problem.Permeability(cf.Subdomain(), c());
const Point &Qf_c = felem.QPoint(q);
int q1 = findQPointID(felem, felem_1, Qf_c);
Tensor K_1 = problem.Permeability(cf.Subdomain(), cf());
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
for (int j = 0; j < felem.dimension(); ++j) {
Scalar phi_j = felem.Value(q, j);
Scalar NDphi_j = (K * felem.Derivative(q, j)) * N;
A_c(i, j) += w * (-0.5 * NDphi_i * phi_j * sign
- 0.5 * phi_i * NDphi_j
+ s * phi_i * phi_j);
}
for (int j = 0; j < felem_1.dimension(); ++j) {
Scalar phi_j = felem_1.Value(q1, j);
Scalar NDphi_j = (K_1 * felem_1.Derivative(q1, j)) * N;
A_cf(i, j) += w * (0.5 * NDphi_i * phi_j * sign
- phi_i * 0.5 * NDphi_j
- s * phi_i * phi_j);
A_fc(j, i) += w * (0.5 * NDphi_i * phi_j
- phi_i * 0.5 * NDphi_j * sign
- s * phi_i * phi_j);
}
}
for (int i = 0; i < felem_1.dimension(); ++i) {
Scalar phi_i = felem_1.Value(q1, i);
Scalar NDphi_i = (K_1 * felem_1.Derivative(q1, i)) * N;
for (int j = 0; j < felem_1.dimension(); ++j) {
Scalar phi_j = felem_1.Value(q1, j);
Scalar NDphi_j = (K_1 * felem_1.Derivative(q1, j)) * N;
A_ff(i, j) += w * (0.5 * NDphi_i * phi_j * sign
+ phi_i * 0.5 * NDphi_j
+ s * phi_i * phi_j);
}
}
}
}
}
}
void DGEllipticAssemble::setExactSolution(Vector &u) const {
u = 0;
Point z[MaxNodalPoints];
for (cell c = u.cells(); c != u.cells_end(); ++c) {
row r = u.find_row(c());
DGElement Elem(disc, u, c);
Elem.NodalPoints(z);
for (int j = 0; j < Elem.dimension(); ++j)
u(r, j) = problem.Solution(z[j]);
}
Accumulate(u);
}
void DGEllipticAssemble::AssembleTransfer(TransferMatrix &TM) const {
TM = 0;
const matrixgraph &cg = TM.CoarseMatrixGraph();
const matrixgraph &fg = TM.FineMatrixGraph();
const Mesh &mesh = fg.GetMesh();
for (cell c = cg.cells(); c != cg.cells_end(); ++c)
for (int k = 0; k < c.Children(); ++k)
TM(c(), c.Child(k))[0] = 1;
}
void DGEllipticAssemble::plotU(Plot &plot, const Vector &u, Meshes &meshes) {
Vector plot_u(u);
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(disc, u, c);
double U = 0;
double area = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
U += w * elem.Value(q, u);
area += w;
}
U *= (1 / area);
plot_u(c(), 0) = U;
}
plot.celldata(plot_u);
plot.vtk_celldata("u");
}
void DGEllipticAssemble::plotUex(Plot &plot, const Vector &u, Meshes &meshes) {
if (!problem.hasExactSolution()) return;
Vector u_ex(u);
setExactSolution(u_ex);
Vector plot_u(u_ex);
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(disc, u, c);
double U = 0;
double area = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
U += w * elem.Value(q, u_ex);
area += w;
}
U *= (1 / area);
plot_u(c(), 0) = U;
}
plot.celldata(plot_u);
plot.vtk_celldata("u_ex");
}
mlmc/src/assemble/DGEllipticAssemble.h 0 → 100644
View file @ 347756ff
#ifndef _DGLAPLACE_H_
#define _DGLAPLACE_H_
#include "EllipticAssemble.h"
#include "dg/DGDiscretization.h"
#include "dg/DGUtils.h"
class DGEllipticAssemble : public EllipticAssemble, public DGUtils {
public:
const DGDiscretization &disc;
const ScalarProblem &problem;
double penalty = 1.0;
int sign = 1;
DGEllipticAssemble(const DGDiscretization &disc, const ScalarProblem &problem);
DGEllipticAssemble(const DGDiscretization &disc, const ScalarProblem &problem,
double penalty, int sign);
Element *getElement(const cell &c, const Vector &u) const override;
FaceElement *getFaceElement(const cell &c, const Vector &u, int face) const override;
Scalar getValue(const cell &c, const Vector &u, Element *elem, int q) const override;
Scalar getFaceValue(const cell &c, const Vector &u,
FaceElement *faceElem, int q, int face) const override;
VectorField getFlux(const cell &c, const Vector &u,
Element *elem, int q) const override;
VectorField getFaceFlux(const cell &c, const Vector &u,
FaceElement *faceElem, int q, int face) const override;
VectorField getGradient(const cell &c, const Vector &u,
Element *elem, int q) const override;
const char *Name() const override;
void Dirichlet(const cell &c, Vector &u) const override;
void Residual(const cell &c, const Vector &u, Vector &r) const override;
void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
void setExactSolution(Vector &u) const override;
void AssembleTransfer(TransferMatrix &TM) const override;
void plotU(Plot &plot, const Vector &u, Meshes &meshes) override;
void plotUex(Plot &plot, const Vector &u, Meshes &meshes) override;
};
#endif
mlmc/src/assemble/DGReactionAssemble.h 0 → 100644
View file @ 347756ff
#ifndef _DG_REACTION_H_
#define _DG_REACTION_H_
#include "dg/DGTAssemble.h"
class DGReactionAssemble : public DGTAssemble {
public:
const ScalarProblem &problem;
double penalty = 1.0;
int sign = 1;
DGReactionAssemble(const DGDiscretization &disc, const ScalarProblem &problem,
Plot &plot) : DGTAssemble(disc, plot), problem(problem) {
ReadConfig(Settings, "penalty", penalty);
ReadConfig(Settings, "sign", sign);
}
const char *Name() const override { return "DGReactionAssemble"; }
void PrintInfo(const Vector &u) override {
Assemble::PrintInfo(u);
Vout (1) << " Problem: " << problem.Name() << endl
<< " Discretization: " << disc.DiscName() << endl
<< endl;
}
double Residual(double t, const Vector &u, Vector &r) const override {
r = 0;
const Vector &u_old = U_old();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
int sd = c.Subdomain();
RowBndValues r_c(r, c);
DGElement elem(disc, u, c);
double kappa = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem.Diffusion(sd, z);
kappa = max(kappa, norm(K));
Scalar f = problem.Load(sd, z);
Scalar U = elem.Value(q, u);
VectorField DU = elem.Derivative(q, u);
VectorField K_DU = K * DU;
Scalar Uold = elem.Value(q, u_old);
Scalar R = problem.Reaction(z, U);
VectorField B = problem.cellConvection(sd, c, z);
for (int i = 0; i < elem.dimension(); ++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);
}
}
for (int f = 0; f < c.Faces(); ++f) {
Scalar s = penalty * c.FaceArea(f);
DGFaceElement felem(disc, u, c, f);
if (checkIfBND(r, c, f)) {
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar BN = problem.faceConvection(sd, c, f, N, z);
Tensor K = problem.Diffusion(c.Subdomain(), z);
Scalar U_diff = felem.Value(q, u)
- problem.Concentration(t, z);
VectorField K_DU = K * felem.Derivative(q, u);
if (BN < 0) {
Scalar F = K_DU * N;
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i =
(K * felem.Derivative(q, i)) * N;
r(c(), i) -= w * (F * phi_i
+ BN * U_diff * phi_i
+ U_diff *
(NDphi_i - s * phi_i));
}
} else {
Scalar F = problem.Flux(2, z) * N;
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
r(c(), i) -= w * F * phi_i;
}
}
}
} else {
cell cf = findNBCell(r, c, f);
if (cf() < c()) continue;
int f1 = findNBFaceID(r, c.Face(f), cf);
DGFaceElement felem_1(disc, r, cf, f1);
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar BN = problem.faceConvection(sd, c, f, N, z);
Tensor K = problem.Diffusion(c.Subdomain(), z);
Scalar U = felem.Value(q, u);
VectorField K_DU = K * felem.Derivative(q, u);
Point Qf_c = felem.QPoint(q);
int q1 = findQPointID(felem, felem_1, Qf_c);
Tensor K_1 = problem.Diffusion(cf.Subdomain(), z);
Scalar U_1 = felem_1.Value(q1, u);
VectorField K_DU_1 = K_1 * felem_1.Derivative(q1, u);
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
r(c(), i) += w * ((s * U - 0.5 * (K_DU * N)) * phi_i
- 0.5 * U * NDphi_i);
r(c(), i) +=
w * ((-s * U_1 - 0.5 * (K_DU_1 * N)) * phi_i
+ 0.5 * U_1 * NDphi_i);
if (BN < 0) {
r(c(), i) -= w * BN * U * phi_i;
r(c(), i) += w * BN * U_1 * phi_i;
}
}
for (int j = 0; j < felem_1.dimension(); ++j) {
Scalar phi_j = felem_1.Value(q1, j);
Scalar NDphi_j = (K_1 * felem_1.Derivative(q1, j)) * N;
r(cf(), j) +=
w * ((0.5 * K_DU_1 * N + s * U_1) * phi_j
+ 0.5 * U_1 * NDphi_j);
r(cf(), j) += w * ((0.5 * K_DU * N - s * U) * phi_j
- 0.5 * U * NDphi_j);
if (BN > 0) {
r(cf(), j) += w * BN * U_1 * phi_j;
r(cf(), j) -= w * BN * U * phi_j;
}
}
}
}
}
}
Collect(r);
return r.norm();
}
void Jacobi(double t, const Vector &u, Matrix &A) const override {
A = 0;
for (cell c = A.cells(); c != A.cells_end(); ++c) {
int sd = c.Subdomain();
DGElement elem(disc, A, c);
DGRowEntries A_c(A, c, c);
double kappa = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem.Diffusion(sd, z);
kappa = max(kappa, norm(K));
Scalar U = elem.Value(q, u);
Scalar DR = problem.D_Reaction(z, U);
VectorField B = problem.cellConvection(sd, c, z);
for (int i = 0; i < elem.dimension(); ++i) {
Scalar U_i = elem.Value(q, i);
VectorField DU_i = elem.Derivative(q, i);
VectorField K_DU_i = K * DU_i;
for (int j = 0; j < elem.dimension(); ++j) {
Scalar U_j = elem.Value(q, j);
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
- DR * U_i * U_j);
}
}
}
BFParts bnd(u.GetMesh(), c);
for (int f = 0; f < c.Faces(); ++f) {
Scalar s = penalty * c.FaceArea(f);
DGFaceElement felem(disc, u, c, f);
if (checkIfBND(u, c, f)) {
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar BN = problem.faceConvection(sd, c, f, N, z);
if (BN < 0) {
Tensor K = problem.Diffusion(c.Subdomain(), z);
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
for (int j = 0; j < felem.dimension(); ++j) {
Scalar phi_j = felem.Value(q, j);
Scalar NDphi_j = (K * felem.Derivative(q, j)) * N;
A_c(i, j) -= w * (NDphi_i * phi_j
+ BN * phi_j * phi_i
+ phi_i *
(NDphi_j - s * phi_j));
}
}
}
}
} else {
cell cf = findNBCell(u, c, f);
if (cf() < c()) continue;
DGRowEntries A_cf(A, c, cf);
DGRowEntries A_fc(A, cf, c);
DGRowEntries A_ff(A, cf, cf);