Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Breaking] Linear Solver Updates #1174

Open
wants to merge 72 commits into
base: develop
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
72 commits
Select commit Hold shift + click to select a range
a1869ad
Initial crack at CG
lroberts36 Sep 5, 2024
fe6a903
Messing around with CG, not working well
lroberts36 Sep 12, 2024
507cbd3
Add CG option to poisson
lroberts36 Sep 12, 2024
38543b0
explicitly set beta = 0 on first iter
lroberts36 Sep 18, 2024
162a0a6
Messing around with different prolongation operators
lroberts36 Sep 24, 2024
216f370
Fix bugs
lroberts36 Sep 24, 2024
83226a3
Merge branch 'develop' into lroberts36/add-cg-solver
lroberts36 Sep 24, 2024
3c81e68
format and lint
lroberts36 Sep 24, 2024
a1bd975
allow for index flattening
lroberts36 Sep 24, 2024
627efe7
Add diagonal preconditioning
lroberts36 Sep 24, 2024
8140938
Include boundary info in sparse pack
lroberts36 Sep 24, 2024
879dcb5
small
lroberts36 Sep 24, 2024
63e58e0
Only smooth on active blocks
lroberts36 Sep 24, 2024
caa959a
Allow switching boundary prolongation method
lroberts36 Sep 24, 2024
0da700f
explicitly impose boundary conditions on fluxes
lroberts36 Sep 24, 2024
f8413a4
small
lroberts36 Sep 24, 2024
843a264
make defaults consistent with older versions
lroberts36 Sep 25, 2024
fbc07fe
Merge branch 'develop' into lroberts36/add-cg-solver
lroberts36 Sep 25, 2024
297d9cc
More freedom in setting prolongation
lroberts36 Sep 25, 2024
4fd9883
More selectability
lroberts36 Sep 25, 2024
bca3a1b
remove MG prolongation operations since they are now user defined
lroberts36 Sep 25, 2024
d554036
format
lroberts36 Sep 25, 2024
1632222
update script
lroberts36 Sep 25, 2024
23aeba3
fix strange compiler error
lroberts36 Sep 26, 2024
cdaee59
changelog
lroberts36 Sep 26, 2024
bd0e875
Merge branch 'develop' into lroberts36/add-cg-solver
lroberts36 Sep 26, 2024
e477e2a
small
lroberts36 Sep 26, 2024
9365bb2
shits compiling
lroberts36 Oct 2, 2024
6fb993b
use solver base class
lroberts36 Oct 2, 2024
4906142
Apparently working stage based setup
lroberts36 Oct 2, 2024
62414a0
switch to type lists
lroberts36 Oct 2, 2024
73ad5bc
start, no where near compiling
lroberts36 Oct 2, 2024
6ebb671
seemingly working staged MG
lroberts36 Oct 3, 2024
dee8241
start on bicg stages, not compiling
lroberts36 Oct 3, 2024
3f9712f
seemingly working staged bicgstab solver
lroberts36 Oct 3, 2024
aa7fe69
format and lint
lroberts36 Oct 3, 2024
548dc7f
Allow user defined prolongation in stage based solvers
lroberts36 Oct 3, 2024
c090262
small
lroberts36 Oct 3, 2024
8dd7962
make base class separate file
lroberts36 Oct 3, 2024
dcd7bef
format
lroberts36 Oct 3, 2024
c678449
moving toward separate def of interior prolongation operators
lroberts36 Oct 8, 2024
75ef64c
refactor user definable block interior prolongation
lroberts36 Oct 8, 2024
6b10143
Remove some more stuff
lroberts36 Oct 8, 2024
96ae149
format
lroberts36 Oct 8, 2024
01917a2
possibly fix non-gcc compilation errors
lroberts36 Oct 8, 2024
c176abb
update docs
lroberts36 Oct 8, 2024
9c7276d
Merge branch 'develop' into lroberts36/add-cg-solver
lroberts36 Oct 8, 2024
a31f34f
small change to doc
lroberts36 Oct 9, 2024
9b0a9ac
Only include solver fields in md_u
lroberts36 Oct 10, 2024
93e08b1
explicitly only include the solution fields
lroberts36 Oct 10, 2024
068abc5
Merge branch 'lroberts36/add-more-careful-container-checking' into lr…
lroberts36 Oct 10, 2024
5f455c9
Merge branch 'lroberts36/add-more-careful-container-checking' into lr…
lroberts36 Oct 10, 2024
6fb49a4
format
lroberts36 Oct 10, 2024
c9619c0
Merge branch 'develop' into lroberts36/add-cg-solver
lroberts36 Oct 29, 2024
456be28
Merge branch 'develop' into lroberts36/add-cg-solver
lroberts36 Nov 19, 2024
5b737e9
add unique ids to staged solvers
lroberts36 Nov 20, 2024
83151d2
small
lroberts36 Nov 20, 2024
5f866e7
clean up macros
lroberts36 Nov 21, 2024
bb49af8
add staged GMG for testing
lroberts36 Nov 21, 2024
d00c887
remove type based solvers
lroberts36 Nov 21, 2024
c882ab2
remove more stuff
lroberts36 Nov 21, 2024
ca1f1f5
remove last thing
lroberts36 Nov 21, 2024
37500d2
closer to resolving all stages
lroberts36 Nov 21, 2024
4789e0e
finish cleanup
lroberts36 Nov 21, 2024
6d1a882
possibly fix compilation
lroberts36 Nov 21, 2024
2478da4
put same MeshData utilities in their own namespace
lroberts36 Nov 25, 2024
4cfc2cc
Merge branch 'develop' into lroberts36/add-cg-solver
lroberts36 Nov 25, 2024
d4a5a26
fix changelog
lroberts36 Nov 25, 2024
0d49d78
Update src/solvers/cg_solver.hpp
lroberts36 Nov 27, 2024
95ed5b2
act on Jonah comments
lroberts36 Nov 27, 2024
d71062e
format
lroberts36 Nov 27, 2024
b5b68c2
Merge branch 'develop' into lroberts36/add-cg-solver
Yurlungur Dec 20, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
More selectability
lroberts36 committed Sep 25, 2024
commit 4fd9883bfbe26cf115f29d2a4911a7a435a0a5e3
6 changes: 2 additions & 4 deletions example/poisson_gmg/poisson_driver.cpp
Original file line number Diff line number Diff line change
@@ -75,7 +75,6 @@ TaskCollection PoissonDriver::MakeTaskCollection(BlockList_t &blocks) {

auto pkg = pmesh->packages.Get("poisson_package");
auto solver = pkg->Param<std::string>("solver");
auto flux_correct = pkg->Param<bool>("flux_correct");
auto use_exact_rhs = pkg->Param<bool>("use_exact_rhs");
auto *mg_solver =
pkg->MutableParam<parthenon::solvers::MGSolver<u, rhs, PoissonEquation>>(
@@ -100,9 +99,8 @@ TaskCollection PoissonDriver::MakeTaskCollection(BlockList_t &blocks) {
if (use_exact_rhs) {
auto copy_exact = tl.AddTask(get_rhs, TF(solvers::utils::CopyData<exact, u>), md);
auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(copy_exact, tl, md, true);
PoissonEquation eqs;
eqs.do_flux_cor = flux_correct;
get_rhs = eqs.Ax<u, rhs>(tl, comm, md);
auto *eqs = pkg->MutableParam<PoissonEquation>("poisson_equation");
get_rhs = eqs->Ax<u, rhs>(tl, comm, md);
}

// Set initial solution guess to zero
101 changes: 91 additions & 10 deletions example/poisson_gmg/poisson_equation.hpp
Original file line number Diff line number Diff line change
@@ -37,6 +37,24 @@ 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;

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");
auto pro_int = pin->GetOrAddString(label, "interior_prolongation", "Linear");
if (pro_int == "Constant") {
prolongation_type = ProlongationType::Constant;
} else if (pro_int == "Linear") {
prolongation_type = ProlongationType::Linear;
} else if (pro_int == "Kwak") {
prolongation_type = ProlongationType::Kwak;
} else {
PARTHENON_FAIL("Invalid user prolongation type.");
}
}

// Add tasks to calculate the result of the matrix A (which is implicitly defined by
// this class) being applied to x_t and store it in field out_t
@@ -166,8 +184,31 @@ 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);
}
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
if (d == 1) return (2.0 + !up_bound) / 4.0;
if (d == -1) return (2.0 + !lo_bound) / 4.0;
if (d == 3) return !up_bound / 4.0;
if (d == -3) return !lo_bound / 4.0;
return 0.0;
}

template <ProlongationType prolongation_type, class... var_ts>
static parthenon::TaskStatus
Prolongate(std::shared_ptr<parthenon::MeshData<Real>> &md) {
ProlongateImpl(std::shared_ptr<parthenon::MeshData<Real>> &md) {
using namespace parthenon;
const int ndim = md->GetMeshPointer()->ndim;
IndexRange ib = md->GetBoundsI(IndexDomain::interior);
@@ -188,7 +229,7 @@ class PoissonEquation {
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),
@@ -197,14 +238,54 @@ class PoissonEquation {
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;
pack(b, n, fk, fj, fi) = pack_coarse(b, n, ck, cj, ci);
//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 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)};

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 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);

}
}
}
} else if constexpr (ProlongationType::Kwak == prolongation_type) {
pack(b, n, fk, fj, fi) = 0.0;
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);
}
}
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);
}
}
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) /= 2.0 * ndim;
}
});
return TaskStatus::complete;
}
16 changes: 3 additions & 13 deletions example/poisson_gmg/poisson_package.cpp
Original file line number Diff line number Diff line change
@@ -78,30 +78,20 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
pkg->UserBoundaryFunctions[BF::outer_x2].push_back(GetBC<X2DIR, BCSide::Outer>());
pkg->UserBoundaryFunctions[BF::outer_x3].push_back(GetBC<X3DIR, BCSide::Outer>());

int max_poisson_iterations = pin->GetOrAddInteger("poisson", "max_iterations", 10000);
pkg->AddParam<>("max_iterations", max_poisson_iterations);

Real diagonal_alpha = pin->GetOrAddReal("poisson", "diagonal_alpha", 0.0);
pkg->AddParam<>("diagonal_alpha", diagonal_alpha);

std::string solver = pin->GetOrAddString("poisson", "solver", "MG");
pkg->AddParam<>("solver", solver);

Real err_tol = pin->GetOrAddReal("poisson", "error_tolerance", 1.e-8);
pkg->AddParam<>("error_tolerance", err_tol);

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

bool flux_correct = pin->GetOrAddBoolean("poisson", "flux_correct", false);
pkg->AddParam<>("flux_correct", flux_correct);
std::string prolong = pin->GetOrAddString("poisson", "boundary_prolongation", "Linear");

std::string prolong = pin->GetOrAddString("poisson", "prolongation", "Linear");
PoissonEquation eq(pin, "poisson");
pkg->AddParam<>("poisson_equation", eq, parthenon::Params::Mutability::Mutable);

PoissonEquation eq;
eq.do_flux_cor = flux_correct;
eq.set_flux_boundary = pin->GetOrAddBoolean("poisson", "set_flux_boundary", false);
eq.include_flux_dx = (prolong == "Constant");
parthenon::solvers::MGParams mg_params(pin, "poisson/solver_params");
parthenon::solvers::MGSolver<u, rhs, PoissonEquation> mg_solver(pkg.get(), mg_params,
eq);
3 changes: 2 additions & 1 deletion src/solvers/bicgstab_solver.hpp
Original file line number Diff line number Diff line change
@@ -151,7 +151,7 @@ class BiCGSTABSolver {
auto copy_rhat0 = tl.AddTask(dependence, TF(CopyData<rhs, rhat0>), md);
auto get_rhat0r_init = DotProduct<rhat0, r>(dependence, tl, &rhat0r, md);
auto get_rhs2 = get_rhat0r_init;
if (params_.relative_residual)
if (params_.relative_residual || params_.print_per_step)
get_rhs2 = DotProduct<rhs, rhs>(dependence, tl, &rhs2, md);
auto initialize = tl.AddTask(
TaskQualifier::once_per_region | TaskQualifier::local_sync,
@@ -172,6 +172,7 @@ class BiCGSTABSolver {
: *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()));
}
return TaskStatus::complete;
},
9 changes: 5 additions & 4 deletions src/solvers/cg_solver.hpp
Original file line number Diff line number Diff line change
@@ -121,7 +121,7 @@ class CGSolver {
auto zero_p = tl.AddTask(dependence, TF(SetToZero<p>), md);
auto copy_r = tl.AddTask(dependence, TF(CopyData<rhs, r>), md);
auto get_rhs2 = none;
if (params_.relative_residual)
if (params_.relative_residual || params_.print_per_step)
get_rhs2 = DotProduct<rhs, rhs>(dependence, tl, &rhs2, md);
auto initialize = tl.AddTask(
TaskQualifier::once_per_region | TaskQualifier::local_sync,
@@ -136,16 +136,17 @@ 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) {
[&](CGSolver *solver, std::shared_ptr<Real> res_tol, bool relative_residual, Mesh *pm) {
Real tol =
relative_residual
? *res_tol * std::sqrt(solver->rhs2.val / pmesh->GetTotalCells())
? *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()));
return TaskStatus::complete;
},
this, params_.residual_tolerance, params_.relative_residual);
this, params_.residual_tolerance, params_.relative_residual, pmesh);
}

// BEGIN ITERATIVE TASKS
2 changes: 1 addition & 1 deletion src/solvers/mg_solver.hpp
Original file line number Diff line number Diff line change
@@ -506,7 +506,7 @@ class MGSolver {
recv_from_coarser, BTF(SetBounds<BoundaryType::gmg_prolongate_recv>), md_comm);
auto prolongate = set_from_coarser;
if (params_.prolongation == "User") {
prolongate = tl.AddTask(set_from_coarser, BTF(&equations::template Prolongate<res_err>), md_comm);
prolongate = eqs_.template Prolongate<res_err>(tl, set_from_coarser, md_comm);
} else {
prolongate = tl.AddTask(set_from_coarser,
BTF(ProlongateBounds<BoundaryType::gmg_prolongate_recv>), md_comm);