Skip to content

Add user specifiable BCs to solvers (+ other small solver updates) #1315

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.

Sign up for GitHub

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

Merged
lroberts36 merged 6 commits into from
Sep 16, 2025
Merged

Add user specifiable BCs to solvers (+ other small solver updates) #1315

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
684a28c
allow for user specifiable BCs, make counters work, add more to base …
lroberts36 Sep 4, 2025
c1a1cf8
format
lroberts36 Sep 4, 2025
2052ad2
Add some strangely missing includes
lroberts36 Sep 4, 2025
b0b45d3
Merge branch 'develop' into lroberts36/add-user-specifiable-BCs-to-so…
lroberts36 Sep 15, 2025
32e4dca
format and lint
lroberts36 Sep 15, 2025
bbb89bb
changelog
lroberts36 Sep 15, 2025
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
59 changes: 32 additions & 27 deletions src/solvers/bicgstab_solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include <utility>
#include <vector>

#include "bvals/comms/bvals_in_one.hpp"
#include "interface/mesh_data.hpp"
#include "interface/meshblock_data.hpp"
#include "interface/state_descriptor.hpp"
Expand All @@ -28,6 +29,7 @@
#include "solvers/solver_utils.hpp"
#include "tasks/tasks.hpp"
#include "utils/type_list.hpp"
#include "utils/reductions.hpp"

namespace parthenon {

Expand Down Expand Up @@ -63,42 +65,34 @@ struct BiCGSTABParams {
}
};

// The equations class must include a template method
struct BiCGSTABSolverCounter {
static inline std::size_t id{0};
};

// The equations_t class must include a template method
//
// template <class x_t, class y_t, class TL_t>
// TaskID Ax(TL_t &tl, TaskID depends_on, std::shared_ptr<MeshData<Real>> &md)
//
// that takes a field associated with x_t and applies
// the matrix A to it and stores the result in y_t.
template <class equations, class preconditioner_t = MGSolver<equations>>
class BiCGSTABSolver : public SolverBase {
using FieldTL = typename equations::IndependentVars;

std::vector<std::string> sol_fields;
// Name of user defined container that should contain information required to
// calculate the matrix part of the matrix vector product
std::string container_base;
// User defined container in which the solution will reside, only needs to contain
// sol_fields
// TODO(LFR): Also allow for an initial guess to come in here
std::string container_u;
// User defined container containing the rhs vector, only needs to contain sol_fields
std::string container_rhs;
template <class equations_t, class preconditioner_t = MGSolver<equations_t>>
class BiCGSTABSolver : public SolverBase, BiCGSTABSolverCounter {
using FieldTL = typename equations_t::IndependentVars;

// Internal containers for solver which create deep copies of sol_fields
std::string container_rhat0, container_v, container_h, container_s;
std::string container_t, container_r, container_p, container_x, container_diag;

static inline std::size_t id{0};
BValOnMDFunc_t BCFunc;

public:
BiCGSTABSolver(const std::string &container_base, const std::string &container_u,
const std::string &container_rhs, ParameterInput *pin,
const std::string &input_block, equations eq_in = equations())
const std::string &input_block, equations_t eq_in = equations_t())
: preconditioner(container_base, container_u, container_rhs, pin, input_block,
eq_in),
container_base(container_base), container_u(container_u),
container_rhs(container_rhs), params_(pin, input_block), iter_counter(0),
eqs_(eq_in) {
SolverBase(container_base, container_u, container_rhs), params_(pin, input_block),
iter_counter(0), eqs_(eq_in) {
FieldTL::IterateTypes(
[this](auto t) { this->sol_fields.push_back(decltype(t)::name()); });
std::string solver_id = "bicgstab" + std::to_string(id++);
Expand All @@ -111,6 +105,17 @@ class BiCGSTABSolver : public SolverBase {
container_p = solver_id + "_p";
container_x = solver_id + "_x";
container_diag = solver_id + "_diag";
if constexpr (has_SetBoundary<equations_t>::value) {
BCFunc = equations_t::SetBoundary;
} else {
BCFunc = ApplyBoundaryConditionsOnCoarseOrFineMD;
}
}

TaskID Ax(TaskList &tl, TaskID dependence, std::shared_ptr<MeshData<Real>> &md_mat,
std::shared_ptr<MeshData<Real>> &md_in,
std::shared_ptr<MeshData<Real>> &md_out) {
return eqs_.Ax(tl, dependence, md_mat, md_in, md_out);
}

TaskID AddSetupTasks(TaskList &tl, TaskID dependence, int partition, Mesh *pmesh) {
Expand All @@ -120,7 +125,7 @@ class BiCGSTABSolver : public SolverBase {
auto partitions = pmesh->GetDefaultBlockPartitions();
auto &md = pmesh->mesh_data.Add(container_base, partitions[partition]);
auto &md_diag = pmesh->mesh_data.Add(container_diag, md, sol_fields);
return tl.AddTask(dependence, &equations::SetDiagonal, &eqs_, md, md_diag);
return tl.AddTask(dependence, &equations_t::SetDiagonal, &eqs_, md, md_diag);
} else {
return dependence;
}
Expand Down Expand Up @@ -217,8 +222,8 @@ class BiCGSTABSolver : public SolverBase {
}

// 2. v <- A u
auto comm =
AddBoundaryExchangeTasks<BoundaryType::any>(precon1, itl, md_u, multilevel);
auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(precon1, itl, md_u,
multilevel, BCFunc);
auto get_v = eqs_.Ax(itl, comm, md_base, md_u, md_v);

// 3. rhat0v <- (rhat0, v)
Expand Down Expand Up @@ -271,8 +276,8 @@ class BiCGSTABSolver : public SolverBase {
}

// 7. t <- A u
auto pre_t_comm =
AddBoundaryExchangeTasks<BoundaryType::any>(precon2, itl, md_u, multilevel);
auto pre_t_comm = AddBoundaryExchangeTasks<BoundaryType::any>(precon2, itl, md_u,
multilevel, BCFunc);
auto get_t = eqs_.Ax(itl, pre_t_comm, md_base, md_u, md_t);

// 8. omega <- (t,s) / (t,t)
Expand Down Expand Up @@ -365,7 +370,7 @@ class BiCGSTABSolver : public SolverBase {
int iter_counter;
AllReduce<Real> rtr, pAp, rhat0v, rhat0r, ts, tt, residual, rhs2;
Real rhat0r_old;
equations eqs_;
equations_t eqs_;
std::string container_;
};

Expand Down
60 changes: 33 additions & 27 deletions src/solvers/cg_solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,11 @@
#include <limits>
#include <memory>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>

#include "bvals/comms/bvals_in_one.hpp"
#include "interface/mesh_data.hpp"
#include "interface/meshblock_data.hpp"
#include "interface/state_descriptor.hpp"
Expand All @@ -29,6 +31,7 @@
#include "solvers/solver_utils.hpp"
#include "tasks/tasks.hpp"
#include "utils/type_list.hpp"
#include "utils/reductions.hpp"

namespace parthenon {

Expand All @@ -54,48 +57,51 @@ struct CGParams {
}
};

// The equations class must include a template method
struct CGSolverCounter {
static inline std::size_t id{0};
};

// The equations_t class must include a template method
//
// template <class x_t, class y_t, class TL_t>
// TaskID Ax(TL_t &tl, TaskID depends_on, std::shared_ptr<MeshData<Real>> &md)
//
// that takes a field associated with x_t and applies
// the matrix A to it and stores the result in y_t.
template <class equations, class preconditioner_t = MGSolver<equations>>
class CGSolver : public SolverBase {
using FieldTL = typename equations::IndependentVars;

std::vector<std::string> sol_fields;
// Name of user defined container that should contain information required to
// calculate the matrix part of the matrix vector product
std::string container_base;
// User defined container in which the solution will reside, only needs to contain
// sol_fields
// TODO(LFR): Also allow for an initial guess to come in here
std::string container_u;
// User defined container containing the rhs vector, only needs to contain sol_fields
std::string container_rhs;
template <class equations_t, class preconditioner_t = MGSolver<equations_t>>
class CGSolver : public SolverBase, CGSolverCounter {
using FieldTL = typename equations_t::IndependentVars;

// Internal containers for solver which create deep copies of sol_fields
std::string container_x, container_r, container_v, container_p;

static inline std::size_t id{0};
BValOnMDFunc_t BCFunc;

public:
CGSolver(const std::string &container_base, const std::string &container_u,
const std::string &container_rhs, ParameterInput *pin,
const std::string &input_block, const equations &eq_in = equations())
: preconditioner(container_base, container_u, container_rhs, pin, input_block,
const std::string &input_block, const equations_t &eq_in = equations_t())
: SolverBase(container_base, container_u, container_rhs),
preconditioner(container_base, container_u, container_rhs, pin, input_block,
eq_in),
container_base(container_base), container_u(container_u),
container_rhs(container_rhs), params_(pin, input_block), iter_counter(0),
eqs_(eq_in) {
params_(pin, input_block), iter_counter(0), eqs_(eq_in) {
FieldTL::IterateTypes(
[this](auto t) { this->sol_fields.push_back(decltype(t)::name()); });
std::string solver_id = "cg" + std::to_string(id++);
container_x = solver_id + "_x";
container_r = solver_id + "_r";
container_v = solver_id + "_v";
container_p = solver_id + "_p";
if constexpr (has_SetBoundary<equations_t>::value) {
BCFunc = equations_t::SetBoundary;
} else {
BCFunc = ApplyBoundaryConditionsOnCoarseOrFineMD;
}
}

TaskID Ax(TaskList &tl, TaskID dependence, std::shared_ptr<MeshData<Real>> &md_mat,
std::shared_ptr<MeshData<Real>> &md_in,
std::shared_ptr<MeshData<Real>> &md_out) {
return eqs_.Ax(tl, dependence, md_mat, md_in, md_out);
}

TaskID AddSetupTasks(TaskList &tl, TaskID dependence, int partition, Mesh *pmesh) {
Expand Down Expand Up @@ -154,7 +160,7 @@ class CGSolver : public SolverBase {
: *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()));
printf("\t0 %e\n", std::sqrt(solver->rhs2.val / pm->GetTotalCells()));
return TaskStatus::complete;
},
this, params_.residual_tolerance, params_.relative_residual, pmesh);
Expand Down Expand Up @@ -199,8 +205,8 @@ class CGSolver : public SolverBase {
this, md_u, md_p);

// 4. v <- A p
auto comm =
AddBoundaryExchangeTasks<BoundaryType::any>(correct_p, itl, md_p, multilevel);
auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(correct_p, itl, md_p,
multilevel, BCFunc);
auto get_v = eqs_.Ax(itl, comm, md_base, md_p, md_v);

// 5. alpha <- r dot u / p dot v (calculate denominator)
Expand Down Expand Up @@ -234,7 +240,7 @@ class CGSolver : public SolverBase {
[&](CGSolver *solver, Mesh *pmesh) {
Real rms_res = std::sqrt(solver->residual.val / pmesh->GetTotalCells());
if (Globals::my_rank == 0 && solver->params_.print_per_step)
printf("\t%i %e\n", solver->iter_counter, rms_res);
printf("\t%i %e\n", solver->iter_counter + 1, rms_res);
return TaskStatus::complete;
},
this, pmesh);
Expand Down Expand Up @@ -273,7 +279,7 @@ class CGSolver : public SolverBase {
int iter_counter;
AllReduce<Real> ru, pAp, residual, rhs2;
Real ru_old;
equations eqs_;
equations_t eqs_;
};

} // namespace solvers
Expand Down
50 changes: 27 additions & 23 deletions src/solvers/mg_solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include <utility>
#include <vector>

#include "bvals/comms/bvals_in_one.hpp"
#include "interface/mesh_data.hpp"
#include "interface/meshblock_data.hpp"
#include "interface/state_descriptor.hpp"
Expand All @@ -31,6 +32,7 @@
#include "tasks/tasks.hpp"
#include "utils/robust.hpp"
#include "utils/type_list.hpp"
#include "utils/reductions.hpp"

namespace parthenon {

Expand Down Expand Up @@ -60,6 +62,9 @@ struct MGParams {
}
};

struct MGSolverCounter {
static inline std::size_t id{0};
};
// The equations_t class must include a template method
//
// template <class x_t, class y_t, class TL_t>
Expand All @@ -75,25 +80,13 @@ struct MGParams {
// That stores the (possibly approximate) diagonal of matrix A in the field
// associated with the type diag_t. This is used for Jacobi iteration.
template <class equations_t, class prolongator_t = ProlongationBlockInteriorDefault>
class MGSolver : public SolverBase {
static inline std::size_t id{0};

class MGSolver : public SolverBase, MGSolverCounter {
public:
using FieldTL = typename equations_t::IndependentVars;

std::vector<std::string> sol_fields;

// Name of user defined container that should contain information required to
// calculate the matrix part of the matrix vector product
std::string container_base;
// User defined container in which the solution will reside, only needs to contain
// sol_fields
// TODO(LFR): Also allow for an initial guess to come in here
std::string container_u;
// User defined container containing the rhs vector, only needs to contain sol_fields
std::string container_rhs;
// Internal containers for solver which create deep copies of sol_fields
std::string container_res_err, container_temp, container_u0, container_diag;
BValOnMDFunc_t BCFunc;

MGSolver(const std::string &container_base, const std::string &container_u,
const std::string &container_rhs, ParameterInput *pin,
Expand All @@ -104,16 +97,26 @@ class MGSolver : public SolverBase {
MGSolver(const std::string &container_base, const std::string &container_u,
const std::string &container_rhs, MGParams params_in,
equations_t eq_in = equations_t(), prolongator_t prol_in = prolongator_t())
: container_base(container_base), container_u(container_u),
container_rhs(container_rhs), params_(params_in), iter_counter(0), eqs_(eq_in),
prolongator_(prol_in) {
: SolverBase(container_base, container_u, container_rhs), params_(params_in),
iter_counter(0), eqs_(eq_in), prolongator_(prol_in) {
FieldTL::IterateTypes(
[this](auto t) { this->sol_fields.push_back(decltype(t)::name()); });
std::string solver_id = "mg" + std::to_string(id++);
container_res_err = solver_id + "_res_err";
container_temp = solver_id + "_temp";
container_u0 = solver_id + "_u0";
container_diag = solver_id + "_diag";
if constexpr (has_SetBoundary<equations_t>::value) {
BCFunc = equations_t::SetBoundary;
} else {
BCFunc = ApplyBoundaryConditionsOnCoarseOrFineMD;
}
}

TaskID Ax(TaskList &tl, TaskID dependence, std::shared_ptr<MeshData<Real>> &md_mat,
std::shared_ptr<MeshData<Real>> &md_in,
std::shared_ptr<MeshData<Real>> &md_out) {
return eqs_.Ax(tl, dependence, md_mat, md_in, md_out);
}

TaskID AddTasks(TaskList &tl, TaskID dependence, const int partition, Mesh *pmesh) {
Expand All @@ -140,7 +143,8 @@ class MGSolver : public SolverBase {
auto &md_res_err = pmesh->mesh_data.Add(container_res_err, md, sol_fields);
auto &md_rhs = pmesh->mesh_data.Add(container_rhs, md, sol_fields);
auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(mg_finest, itl, md_u,
pmesh->multilevel);
pmesh->multilevel, BCFunc);

auto calc_pointwise_res = eqs_.Ax(itl, comm, md, md_u, md_res_err);
calc_pointwise_res =
itl.AddTask(calc_pointwise_res, TF(AddFieldsAndStoreInteriorSelect<FieldTL>),
Expand Down Expand Up @@ -288,8 +292,8 @@ class MGSolver : public SolverBase {
auto &md_rhs = pmesh->mesh_data.Add(container_rhs, partitions[partition], sol_fields);
auto &md_diag = pmesh->mesh_data.Add(container_diag, md_base, sol_fields);

auto comm =
AddBoundaryExchangeTasks<comm_boundary>(depends_on, tl, md_in, multilevel);
auto comm = AddBoundaryExchangeTasks<comm_boundary>(depends_on, tl, md_in, multilevel,
BCFunc);
auto mat_mult = eqs_.Ax(tl, comm, md_base, md_in, md_out);
return tl.AddTask(mat_mult, TF(&MGSolver::Jacobi), this, md_rhs, md_out, md_diag,
md_in, md_out, omega);
Expand Down Expand Up @@ -436,7 +440,7 @@ class MGSolver : public SolverBase {
// calling Ax. That being said, at least in one case commenting this line out
// didn't seem to impact the solution.
set_from_finer = AddBoundaryExchangeTasks<BoundaryType::gmg_same>(
set_from_finer, tl, md_u, multilevel);
set_from_finer, tl, md_u, multilevel, BCFunc);
set_from_finer =
tl.AddTask(set_from_finer, BTF(CopyData<FieldTL, true>), md_u, md_u0);
// This should set the rhs only in blocks that correspond to interior nodes, the
Expand All @@ -460,7 +464,7 @@ class MGSolver : public SolverBase {
if (level > min_level) {
// 3. Communicate same level boundaries so that u is up to date everywhere
auto comm_u = AddBoundaryExchangeTasks<BoundaryType::gmg_same>(pre_smooth, tl, md_u,
multilevel);
multilevel, BCFunc);

// 4. Caclulate residual and store in communication field
auto residual = eqs_.Ax(tl, comm_u, md, md_u, md_temp);
Expand Down Expand Up @@ -515,7 +519,7 @@ class MGSolver : public SolverBase {
// This is required to make sure boundaries of res_err are up to date before
// prolongation
auto boundary = AddBoundaryExchangeTasks<BoundaryType::gmg_same>(
copy_over, tl, md_res_err, multilevel);
copy_over, tl, md_res_err, multilevel, BCFunc);
last_task = tl.AddTask(
boundary, BTF(SendBoundBufs<BoundaryType::gmg_prolongate_send>), md_res_err);
}
Expand Down
Loading
Loading