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

deploy commsplit in assemble

parent b117b801
......@@ -164,7 +164,7 @@ double DGLinearTransportAssemble::Energy(const Vector &u) const {
energy += w * (U * U);
}
}
return 0.5 * PPM->Sum(energy);
return 0.5 * PPM->SumOnCommSplit(energy, u.GetMesh().CommSplit());
}
void DGNonLinearTransportAssemble::Energy(const cell &c,
......@@ -189,7 +189,7 @@ double DGLinearTransportAssemble::Error(double t, const Vector &u) const {
err += w * (U - Sol) * (U - Sol);
}
}
return sqrt(PPM->Sum(err));
return sqrt(PPM->SumOnCommSplit(err, u.GetMesh().CommSplit()));
}
std::pair<double, double> DGLinearTransportAssemble::InflowOutflow(const Vector &u) const {
......@@ -219,8 +219,8 @@ std::pair<double, double> DGLinearTransportAssemble::InflowOutflow(const Vector
}
}
}
inflow = PPM->Sum(inflow);
outflow = PPM->Sum(outflow);
inflow = PPM->SumOnCommSplit(inflow, u.GetMesh().CommSplit());
outflow = PPM->SumOnCommSplit(outflow, u.GetMesh().CommSplit());
if (inflow < 1e-7)
inflow = 0.0;
if (outflow < 1e-7)
......@@ -249,7 +249,7 @@ double DGLinearTransportAssemble::Mass(const Vector &u) const {
e += w * U;
}
}
return PPM->Sum(e);
return PPM->SumOnCommSplit(e, u.GetMesh().CommSplit());
}
void DGLinearTransportAssemble::SetInitialValue(Vector &u) const {
......
......@@ -37,7 +37,7 @@ double DGEllipticAssemble::Energy(const Vector &u) const {
energy += w * K_DU * DU;
}
}
return sqrt(PPM->Sum(energy));
return sqrt(PPM->SumOnCommSplit(energy, u.GetMesh().CommSplit()));
}
void DGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) const {
......@@ -228,7 +228,7 @@ double DGEllipticAssemble::EnergyError(const Vector &u) const {
err += w * (IK * diff) * diff;
}
}
return sqrt(PPM->Sum(err));
return sqrt(PPM->SumOnCommSplit(err, u.GetMesh().CommSplit()));
}
double DGEllipticAssemble::L2(const Vector &u) const {
......@@ -241,7 +241,7 @@ double DGEllipticAssemble::L2(const Vector &u) const {
l2 += w * U * U;
}
}
return sqrt(PPM->Sum(l2));
return sqrt(PPM->SumOnCommSplit(l2, u.GetMesh().CommSplit()));
}
double DGEllipticAssemble::H1(const Vector &u) const {
......@@ -255,7 +255,7 @@ double DGEllipticAssemble::H1(const Vector &u) const {
err += w * U * U + w * DU * DU;
}
}
return sqrt(PPM->Sum(err));
return sqrt(PPM->SumOnCommSplit(err, u.GetMesh().CommSplit()));
}
double DGEllipticAssemble::L2Error(const Vector &u) const {
......@@ -269,7 +269,7 @@ double DGEllipticAssemble::L2Error(const Vector &u) const {
err += w * (U - Sol) * (U - Sol);
}
}
return sqrt(PPM->Sum(err));
return sqrt(PPM->SumOnCommSplit(err, u.GetMesh().CommSplit()));
}
double DGEllipticAssemble::L2CellAverageError(const Vector &u) const {
......@@ -286,7 +286,7 @@ double DGEllipticAssemble::L2CellAverageError(const Vector &u) const {
}
err += w * (U - Sol) * (U - Sol);
}
return sqrt(PPM->Sum(err));
return sqrt(PPM->SumOnCommSplit(err, u.GetMesh().CommSplit()));
}
double DGEllipticAssemble::MaxError(const Vector &u) const {
......@@ -298,7 +298,7 @@ double DGEllipticAssemble::MaxError(const Vector &u) const {
err = max(err, abs(p - problem->Solution(elem.QPoint(q))));
}
}
return PPM->Max(err);
return PPM->Max(err, u.GetMesh().CommSplit());
}
double DGEllipticAssemble::FluxError(const Vector &u) const {
......@@ -329,7 +329,7 @@ double DGEllipticAssemble::FluxError(const Vector &u) const {
}
}
}
return sqrt(PPM->Sum(flux_error));
return sqrt(PPM->SumOnCommSplit(flux_error, u.GetMesh().CommSplit()));
}
double DGEllipticAssemble::FaceError(const Vector &u) const {
......@@ -345,7 +345,7 @@ double DGEllipticAssemble::FaceError(const Vector &u) const {
}
}
}
return sqrt(PPM->Sum(face_error));
return sqrt(PPM->SumOnCommSplit(face_error, u.GetMesh().CommSplit()));
}
FluxPair DGEllipticAssemble::InflowOutflow(const Vector &u) const {
......@@ -367,7 +367,8 @@ FluxPair DGEllipticAssemble::InflowOutflow(const Vector &u) const {
}
}
}
return {PPM->Sum(inflow), PPM->Sum(outflow)};
return {PPM->SumOnCommSplit(inflow, u.GetMesh().CommSplit()),
PPM->SumOnCommSplit(outflow, u.GetMesh().CommSplit())};
}
FluxPair DGEllipticAssemble::PrescribedInflowOutflow(const Vector &u) const {
......@@ -388,7 +389,8 @@ FluxPair DGEllipticAssemble::PrescribedInflowOutflow(const Vector &u) const {
}
}
}
return {PPM->Sum(inflow), PPM->Sum(outflow)};
return {PPM->SumOnCommSplit(inflow, u.GetMesh().CommSplit()),
PPM->SumOnCommSplit(outflow, u.GetMesh().CommSplit())};
}
FluxPair DGEllipticAssemble::OutflowLeftRight(const Vector &u) const {
......@@ -410,7 +412,8 @@ FluxPair DGEllipticAssemble::OutflowLeftRight(const Vector &u) const {
}
}
}
return {PPM->Sum(outflowLeft), PPM->Sum(outflowRight)};
return {PPM->SumOnCommSplit(outflowLeft, u.GetMesh().CommSplit()),
PPM->SumOnCommSplit(outflowRight, u.GetMesh().CommSplit())};
}
double DGEllipticAssemble::GoalFunctional(const Vector &u) const {
......@@ -433,7 +436,7 @@ double DGEllipticAssemble::GoalFunctional(const Vector &u) const {
}
}
}
return PPM->Sum(goal);
return PPM->SumOnCommSplit(goal, u.GetMesh().CommSplit());
}
void DGEllipticAssemble::SetExactSolution(Vector &u) const {
......
......@@ -282,7 +282,8 @@ std::pair<double, double> DGReactionAssemble::InFlowOutFlowRate(const Vector &u)
}
}
}
return {PPM->Sum(inflow), PPM->Sum(outflow)};
return {PPM->SumOnCommSplit(inflow, u.GetMesh().CommSplit()),
PPM->SumOnCommSplit(outflow, u.GetMesh().CommSplit())};
}
void DGReactionAssemble::SetInitialValue(Vector &u) const {
......@@ -308,7 +309,7 @@ double DGReactionAssemble::Mass(const Vector &u) const {
e += w * U;
}
}
return PPM->Sum(e);
return PPM->SumOnCommSplit(e, u.GetMesh().CommSplit());
}
void DGReactionAssemble::VtkPlotting_cell(double t,
......
#include "PGReactionAssemble.hpp"
#include "Plotting.hpp"
void PGReactionAssemble::Dirichlet(double t, Vector &u) const {
u.ClearDirichletFlags();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
RowBndValues u_c(u, c);
if (!u_c.onBnd()) continue;
int sd = c.Subdomain();
ScalarElement E(*disc, u, c);
for (int i = 0; i < c.Faces(); ++i) {
if (u_c.bc(i) == -1) continue;
ScalarFaceElement FE(*disc, u, c, i);
Scalar F = problem->FaceConvection(c, i, FE.QNormal(0),
FE.QPoint(0));
if (F > -1e-6) continue;
for (int j = 0; j < disc->NodalPointsOnFace(c, i); ++j) {
int k = disc->NodalPointOnFace(c, i, j);
u_c(k) = problem->Concentration(t, E.NodalPoint(k));
u_c.D(k) = true;
}
}
u.ClearDirichletFlags();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
RowBndValues u_c(u, c);
if (!u_c.onBnd()) continue;
int sd = c.Subdomain();
ScalarElement E(*disc, u, c);
for (int i = 0; i < c.Faces(); ++i) {
if (u_c.bc(i) == -1) continue;
ScalarFaceElement FE(*disc, u, c, i);
Scalar F = problem->FaceConvection(c, i, FE.QNormal(0),
FE.QPoint(0));
if (F > -1e-6) continue;
for (int j = 0; j < disc->NodalPointsOnFace(c, i); ++j) {
int k = disc->NodalPointOnFace(c, i, j);
u_c(k) = problem->Concentration(t, E.NodalPoint(k));
u_c.D(k) = true;
}
}
u.DirichletConsistent();
}
u.DirichletConsistent();
}
double PGReactionAssemble::Residual(double t, 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(*disc, u, c);
RowBndValues r_c(r, c);
double delta_h = Delta(c);
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);
Tensor K = problem->Diffusion(E.QPoint(q));
Scalar f = problem->Load(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));
}
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(*disc, 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));
if (F < 1e-6) continue;
Scalar U = 0;
for (unsigned int i = 0; i < FE.size(); ++i) {
Scalar U_i = FE.Value(q, i);
U += u(FE[i], 0) * U_i;
}
for (unsigned int j = 0; j < FE.size(); ++j) {
Scalar U_j = FE.Value(q, j);
r_f(j) += w * F * U * U_j;
}
}
}
r = 0;
const Vector &u_old = U_old();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
ScalarElement E(*disc, u, c);
RowBndValues r_c(r, c);
double delta_h = Delta(c);
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);
Tensor K = problem->Diffusion(E.QPoint(q));
Scalar f = problem->Load(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));
}
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(*disc, 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));
if (F < 1e-6) continue;
Scalar U = 0;
for (unsigned int i = 0; i < FE.size(); ++i) {
Scalar U_i = FE.Value(q, i);
U += u(FE[i], 0) * U_i;
}
for (unsigned int j = 0; j < FE.size(); ++j) {
Scalar U_j = FE.Value(q, j);
r_f(j) += w * F * U * U_j;
}
}
}
}
r.ClearDirichletValues();
r.Collect();
return r.norm();
}
r.ClearDirichletValues();
r.Collect();
return r.norm();
}
void PGReactionAssemble::Jacobi(double t, const Vector &u, Matrix &A) const {
A = 0;
for (cell c = A.cells(); c != A.cells_end(); ++c) {
ScalarElement E(*disc, u, c);
RowEntries A_c(A, E);
double delta_h = Delta(c);
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
Scalar U = E.Value(q, u);
Tensor K = problem->Diffusion(E.QPoint(q));
VectorField C = problem->CellConvection(
c, E.QPoint(q));
Scalar DR = problem->DerivativeReaction(E.QPoint(q), U);
for (unsigned int i = 0; i < E.size(); ++i) {
VectorField DU_i = E.Derivative(q, i);
Scalar U_i = E.Value(q, i);
for (unsigned int j = 0; j < E.size(); ++j) {
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));
}
}
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(*disc, u, c, f);
RowEntries A_f(A, 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));
if (F < 1e-6) continue;
for (unsigned int i = 0; i < FE.size(); ++i) {
Scalar U_i = FE.Value(q, i);
for (unsigned int j = 0; j < FE.size(); ++j) {
Scalar U_j = FE.Value(q, j);
A_f(i, j) += w * F * U_i * U_j;
}
}
}
A = 0;
for (cell c = A.cells(); c != A.cells_end(); ++c) {
ScalarElement E(*disc, u, c);
RowEntries A_c(A, E);
double delta_h = Delta(c);
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
Scalar U = E.Value(q, u);
Tensor K = problem->Diffusion(E.QPoint(q));
VectorField C = problem->CellConvection(
c, E.QPoint(q));
Scalar DR = problem->DerivativeReaction(E.QPoint(q), U);
for (unsigned int i = 0; i < E.size(); ++i) {
VectorField DU_i = E.Derivative(q, i);
Scalar U_i = E.Value(q, i);
for (unsigned int j = 0; j < E.size(); ++j) {
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));
}
}
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(*disc, u, c, f);
RowEntries A_f(A, 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));
if (F < 1e-6) continue;
for (unsigned int i = 0; i < FE.size(); ++i) {
Scalar U_i = FE.Value(q, i);
for (unsigned int j = 0; j < FE.size(); ++j) {
Scalar U_j = FE.Value(q, j);
A_f(i, j) += w * F * U_i * U_j;
}
}
}
}
}
A.ClearDirichletValues();
}
A.ClearDirichletValues();
}
std::pair<double, double> PGReactionAssemble::InFlowOutFlowRate(const Vector &u) const {
double inflow = 0;
double outflow = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
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(*disc, u, c, f);
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 U = FE.Value(q, u);
if (F > 0) outflow += w * F * U;
else inflow += w * F * U;
}
}
double inflow = 0;
double outflow = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
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(*disc, u, c, f);
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 U = FE.Value(q, u);
if (F > 0) outflow += w * F * U;
else inflow += w * F * U;
}
}
return std::pair<double, double>(PPM->Sum(inflow), PPM->Sum(outflow));
}
return std::pair<double, double>(PPM->SumOnCommSplit(inflow, u.GetMesh().CommSplit()),
PPM->SumOnCommSplit(outflow, u.GetMesh().CommSplit()));
}
void PGReactionAssemble::SetInitialValue(Vector &u) const {
for (cell c = u.cells(); c != u.cells_end(); ++c) {
ScalarElement E(*disc, u, c);
RowBndValues u_c(u, c);
for (unsigned int i = 0; i < E.size(); ++i)
u_c(i) = problem->Concentration(0, E[i]());
}
for (cell c = u.cells(); c != u.cells_end(); ++c) {
ScalarElement E(*disc, u, c);
RowBndValues u_c(u, c);
for (unsigned int i = 0; i < E.size(); ++i)
u_c(i) = problem->Concentration(0, E[i]());
}
}
double PGReactionAssemble::Mass(const Vector &u) const {
Scalar e = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
ScalarElement E(*disc, u, c);
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
Scalar U = E.Value(q, u);
e += w * U;
}
Scalar e = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
ScalarElement E(*disc, u, c);
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
Scalar U = E.Value(q, u);
e += w * U;
}
return PPM->Sum(e);
}
return PPM->SumOnCommSplit(e, u.GetMesh().CommSplit());
}
void PGReactionAssemble::VtkPlotting_cell(double t, const Vector &u, char *filename) const {
mpp::plot("C") << u << mpp::save_plot(NumberName("C", filename, step));
void PGReactionAssemble::VtkPlotting_cell(double t,
const Vector &u,
char *filename) const {
mpp::plot("C") << u << mpp::save_plot(NumberName("C", filename, step));
}
double PGReactionAssemble::Delta(const cell &c) const {
double h = diam(c);
Tensor K = problem->Diffusion(c());
VectorField C = problem->CellConvection(c, c());
double Pe = 0.5 * h * norm(C) / norm(K);
if (Pe > 1) return h * delta;
return h * h * delta / norm(K);
double h = diam(c);
Tensor K = problem->Diffusion(c());
VectorField C = problem->CellConvection(c, c());
double Pe = 0.5 * h * norm(C) / norm(K);
if (Pe > 1) return h * delta;
return h * h * delta / norm(K);
}
double PGReactionAssemble::diam(const cell &c) const {
double h = 0;
for (int i = 0; i < c.Corners(); ++i)
for (int j = i + 1; j < c.Corners(); ++j)
h = max(h, dist(c[i], c[j]));
return h;
double h = 0;
for (int i = 0; i < c.Corners(); ++i)
for (int j = i + 1; j < c.Corners(); ++j)
h = max(h, dist(c[i], c[j]));
return h;
}
......