Skip to content

Commit

Permalink
format
Browse files Browse the repository at this point in the history
  • Loading branch information
@lroberts36
lroberts36 committed Sep 25, 2024
1 parent bca3a1b commit d554036
Show file tree
Hide file tree
Showing 4 changed files with 88 additions and 63 deletions.
38 changes: 25 additions & 13 deletions example/poisson_gmg/plot_convergence.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,29 +13,41 @@

import numpy as np
import glob
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
import subprocess

plt.style.use('tableau-colorblind10')
plt.style.use("tableau-colorblind10")

solver = "BiCGSTAB"
difco = 1e6
for bound_pro in ["Constant", "Linear"]:
for interior_pro in ["Constant", "OldLinear"]:
for interior_pro in ["Constant", "OldLinear"]:

p = subprocess.run(["./poisson-gmg-example", "-i", "parthinput.poisson",
"poisson/solver=" + solver,
"poisson/interior_D=" + str(difco),
"poisson/prolongation=" + bound_pro,
"poisson/solver_params/prolongation=" + interior_pro], capture_output = True, text = True)
dat = np.genfromtxt(p.stdout.splitlines())
p = subprocess.run(
[
"./poisson-gmg-example",
"-i",
"parthinput.poisson",
"poisson/solver=" + solver,
"poisson/interior_D=" + str(difco),
"poisson/prolongation=" + bound_pro,
"poisson/solver_params/prolongation=" + interior_pro,
],
capture_output=True,
text=True,
)
dat = np.genfromtxt(p.stdout.splitlines())

plt.semilogy(dat[:, 0], dat[:, 1], label=solver + "_" + str(difco) + "_" + bound_pro + "_" + interior_pro)
plt.semilogy(
dat[:, 0],
dat[:, 1],
label=solver + "_" + str(difco) + "_" + bound_pro + "_" + interior_pro,
)


plt.legend(loc = 'upper right')
plt.ylim([1.e-14, 1e2])
plt.legend(loc="upper right")
plt.ylim([1.0e-14, 1e2])
plt.xlim([0, 40])
plt.xlabel("# of V-cycles")
plt.ylabel("RMS Residual")
plt.savefig("convergence_1e6.pdf")
plt.savefig("convergence_1e6.pdf")
101 changes: 57 additions & 44 deletions example/poisson_gmg/poisson_equation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,21 +37,22 @@ class PoissonEquation {
bool do_flux_cor = false;
bool set_flux_boundary = false;
bool include_flux_dx = false;
enum class ProlongationType { Constant, Linear, Kwak };
ProlongationType prolongation_type = ProlongationType::Constant;
enum class ProlongationType { Constant, Linear, Kwak };
ProlongationType prolongation_type = ProlongationType::Constant;

PoissonEquation(parthenon::ParameterInput *pin, const std::string &label) {
do_flux_cor = pin->GetOrAddBoolean(label, "flux_correct", false);
set_flux_boundary = pin->GetOrAddBoolean(label, "set_flux_boundary", false);
include_flux_dx = (pin->GetOrAddString(label, "boundary_prolongation", "Linear") == "Constant");
include_flux_dx =
(pin->GetOrAddString(label, "boundary_prolongation", "Linear") == "Constant");
auto pro_int = pin->GetOrAddString(label, "interior_prolongation", "Linear");
if (pro_int == "Constant") {
if (pro_int == "Constant") {
prolongation_type = ProlongationType::Constant;
} else if (pro_int == "Linear") {
} else if (pro_int == "Linear") {
prolongation_type = ProlongationType::Linear;
} else if (pro_int == "Kwak") {
} else if (pro_int == "Kwak") {
prolongation_type = ProlongationType::Kwak;
} else {
} else {
PARTHENON_FAIL("Invalid user prolongation type.");
}
}
Expand Down Expand Up @@ -186,16 +187,19 @@ class PoissonEquation {
template <class... var_ts>
parthenon::TaskID Prolongate(parthenon::TaskList &tl, parthenon::TaskID depends_on,
std::shared_ptr<parthenon::MeshData<Real>> &md) {
if (prolongation_type == ProlongationType::Constant) {
return tl.AddTask(depends_on, ProlongateImpl<ProlongationType::Constant, var_ts...>, md);
} else if (prolongation_type == ProlongationType::Linear) {
return tl.AddTask(depends_on, ProlongateImpl<ProlongationType::Linear, var_ts...>, md);
} else if (prolongation_type == ProlongationType::Kwak) {
return tl.AddTask(depends_on, ProlongateImpl<ProlongationType::Kwak, var_ts...>, md);
if (prolongation_type == ProlongationType::Constant) {
return tl.AddTask(depends_on, ProlongateImpl<ProlongationType::Constant, var_ts...>,
md);
} else if (prolongation_type == ProlongationType::Linear) {
return tl.AddTask(depends_on, ProlongateImpl<ProlongationType::Linear, var_ts...>,
md);
} else if (prolongation_type == ProlongationType::Kwak) {
return tl.AddTask(depends_on, ProlongateImpl<ProlongationType::Kwak, var_ts...>,
md);
}
return depends_on;
}

KOKKOS_FORCEINLINE_FUNCTION
static Real LinearFactor(int d, bool lo_bound, bool up_bound) {
if (d == 0) return 1.0; // Indicates this dimension is not included
Expand All @@ -206,6 +210,15 @@ class PoissonEquation {
return 0.0;
}

KOKKOS_FORCEINLINE_FUNCTION
static Real QuadraticFactor(int d) {
if (d == 0) return 1.0; // Indicates this dimension is not included
if (d == 1 || d == -1) return 30.0 / 32.0;
if (d == 3 || d == -3) return 5.0 / 32.0;
if (d == 5 || d == -5) return -3.0 / 32.0;
return 0.0;
}

template <ProlongationType prolongation_type, class... var_ts>
static parthenon::TaskStatus
ProlongateImpl(std::shared_ptr<parthenon::MeshData<Real>> &md) {
Expand All @@ -222,66 +235,66 @@ class PoissonEquation {

int nblocks = md->NumBlocks();
std::vector<bool> include_block(nblocks, true);
for (int b = 0; b < nblocks; ++b) {
include_block[b] = md->grid.logical_level == md->GetBlockData(b)->GetBlockPointer()->loc.level();
for (int b = 0; b < nblocks; ++b) {
include_block[b] =
md->grid.logical_level == md->GetBlockData(b)->GetBlockPointer()->loc.level();
}
const auto desc = parthenon::MakePackDescriptor<var_ts...>(md.get());
const auto desc_coarse = parthenon::MakePackDescriptor<var_ts...>(md.get(), {}, {PDOpt::Coarse});
const auto desc_coarse =
parthenon::MakePackDescriptor<var_ts...>(md.get(), {}, {PDOpt::Coarse});
auto pack = desc.GetPack(md.get(), include_block);
auto pack_coarse = desc_coarse.GetPack(md.get(), include_block);

parthenon::par_for(
"Prolongate", 0, pack.GetNBlocks() - 1,
pack.GetLowerBoundHost(0), pack.GetUpperBoundHost(0),
kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
KOKKOS_LAMBDA(const int b, const int n, const int fk, const int fj, const int fi) {
"Prolongate", 0, pack.GetNBlocks() - 1, pack.GetLowerBoundHost(0),
pack.GetUpperBoundHost(0), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
KOKKOS_LAMBDA(const int b, const int n, const int fk, const int fj,
const int fi) {
const int ck = (ndim > 2) ? (fk - kb.s) / 2 + ckb.s : ckb.s;
const int cj = (ndim > 1) ? (fj - jb.s) / 2 + cjb.s : cjb.s;
const int ci = (ndim > 0) ? (fi - ib.s) / 2 + cib.s : cib.s;
const int fok = (fk - kb.s) % 2;
const int foj = (fj - jb.s) % 2;
const int foi = (fi - ib.s) % 2;
const bool bound[6]{
pack.IsPhysicalBoundary(b, 0, 0, -1) && (ib.s == fi),
pack.IsPhysicalBoundary(b, 0, 0, 1) && (ib.e == fi),
pack.IsPhysicalBoundary(b, 0, -1, 0) && (jb.s == fj),
pack.IsPhysicalBoundary(b, 0, 1, 0) && (jb.e == fj),
pack.IsPhysicalBoundary(b, -1, 0, 0) && (kb.s == fk),
pack.IsPhysicalBoundary(b, 1, 0, 0) && (kb.e == fk)};

const bool bound[6]{pack.IsPhysicalBoundary(b, 0, 0, -1) && (ib.s == fi),
pack.IsPhysicalBoundary(b, 0, 0, 1) && (ib.e == fi),
pack.IsPhysicalBoundary(b, 0, -1, 0) && (jb.s == fj),
pack.IsPhysicalBoundary(b, 0, 1, 0) && (jb.e == fj),
pack.IsPhysicalBoundary(b, -1, 0, 0) && (kb.s == fk),
pack.IsPhysicalBoundary(b, 1, 0, 0) && (kb.e == fk)};

if constexpr (ProlongationType::Constant == prolongation_type) {
pack(b, n, fk, fj, fi) = pack_coarse(b, n, ck, cj, ci);
} else if constexpr (ProlongationType::Linear == prolongation_type) {
pack(b, n, fk, fj, fi) = 0.0;
for (int ok = -(ndim > 2); ok < 1 + (ndim > 2); ++ok) {
for (int oj = -(ndim > 1); oj < 1 + (ndim > 1); ++oj) {
for (int oi = -(ndim > 0); oi < 1 + (ndim > 0); ++oi) {
const int dx3 = (ndim > 2) ? 4 * ok - (2 * fok - 1) : 0;
const int dx2 = (ndim > 1) ? 4 * oj - (2 * foj - 1) : 0;
const int dx3 = (ndim > 2) ? 4 * ok - (2 * fok - 1) : 0;
const int dx2 = (ndim > 1) ? 4 * oj - (2 * foj - 1) : 0;
const int dx1 = 4 * oi - (2 * foi - 1);
pack(b, n, fk, fj, fi) += LinearFactor(dx1, bound[0], bound[1])
* LinearFactor(dx2, bound[2], bound[3])
* LinearFactor(dx3, bound[4], bound[5])
* pack_coarse(b, n, ck + ok, cj + oj, ci + oi);

pack(b, n, fk, fj, fi) += LinearFactor(dx1, bound[0], bound[1]) *
LinearFactor(dx2, bound[2], bound[3]) *
LinearFactor(dx3, bound[4], bound[5]) *
pack_coarse(b, n, ck + ok, cj + oj, ci + oi);
}
}
}
} else if constexpr (ProlongationType::Kwak == prolongation_type) {
} else if constexpr (ProlongationType::Kwak == prolongation_type) {
pack(b, n, fk, fj, fi) = 0.0;
if (ndim > 2 && !bound[4 + fok]) {
if (ndim > 2 && !bound[4 + fok]) {
for (int ok = fok - 1; ok <= fok; ++ok) {
pack(b, n, fk, fj, fi) += pack_coarse(b, n, ck + ok, cj, ci);
pack(b, n, fk, fj, fi) += pack_coarse(b, n, ck + ok, cj, ci);
}
}
if (ndim > 1 && !bound[2 + foj]) {
if (ndim > 1 && !bound[2 + foj]) {
for (int oj = foj - 1; oj <= foj; ++oj) {
pack(b, n, fk, fj, fi) += pack_coarse(b, n, ck, cj + oj, ci);
pack(b, n, fk, fj, fi) += pack_coarse(b, n, ck, cj + oj, ci);
}
}
if (ndim > 0 && !bound[foi]) {
if (ndim > 0 && !bound[foi]) {
for (int oi = foi - 1; oi <= foi; ++oi) {
pack(b, n, fk, fj, fi) += pack_coarse(b, n, ck, cj, ci + oi);
pack(b, n, fk, fj, fi) += pack_coarse(b, n, ck, cj, ci + oi);
}
}
pack(b, n, fk, fj, fi) /= 2.0 * ndim;
Expand Down
2 changes: 1 addition & 1 deletion example/poisson_gmg/poisson_package.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {

bool use_exact_rhs = pin->GetOrAddBoolean("poisson", "use_exact_rhs", false);
pkg->AddParam<>("use_exact_rhs", use_exact_rhs);

std::string prolong = pin->GetOrAddString("poisson", "boundary_prolongation", "Linear");

PoissonEquation eq(pin, "poisson");
Expand Down
10 changes: 5 additions & 5 deletions src/solvers/cg_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,11 +136,11 @@ class CGSolver {
if (params_.print_per_step && Globals::my_rank == 0) {
initialize = tl.AddTask(
TaskQualifier::once_per_region, initialize, "print to screen",
[&](CGSolver *solver, std::shared_ptr<Real> res_tol, bool relative_residual, Mesh *pm) {
Real tol =
relative_residual
? *res_tol * std::sqrt(solver->rhs2.val / pm->GetTotalCells())
: *res_tol;
[&](CGSolver *solver, std::shared_ptr<Real> res_tol, bool relative_residual,
Mesh *pm) {
Real tol = relative_residual
? *res_tol * std::sqrt(solver->rhs2.val / pm->GetTotalCells())
: *res_tol;
printf("# [0] v-cycle\n# [1] rms-residual (tol = %e) \n# [2] rms-error\n",
tol);
printf("0 %e\n", std::sqrt(solver->rhs2.val / pm->GetTotalCells()));
Expand Down

0 comments on commit d554036

Please sign in to comment.