Skip to content

Commit

Permalink
Segregated momentum solve functionality (Exawind#623)
Browse files Browse the repository at this point in the history 
* Segregated momentum solve algorithm solves u, v and w momentum equations separately.
  Reference PR Exawind#394: Exawind@7ee3a5a

* Feature can be enabled by setting,
  velocity_diffusion.use_tensor_operator = false, and
  velocity_diffusion.use_segregated_operator = true

* Multicomponent solve can enabled by setting both input parameters to false.

* Supports mesh mapping

* Created a new regression test directory
  Test name: "abl_godunov_segregated_velocity_solve"

** Note: The original tensor (default) and multicomponent solve workflows are not affected by the changes.
  • Loading branch information
@sbidadi9
sbidadi9 authored Jun 20, 2022
1 parent 9d851fc commit 7b68af2
Show file tree
Hide file tree
Showing 4 changed files with 441 additions and 5 deletions.
282 changes: 277 additions & 5 deletions amr-wind/equation_systems/icns/icns_diffusion.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -285,6 +285,254 @@ protected:
std::unique_ptr<amrex::MLABecLaplacian> m_applier_scalar;
};

class ICNSDiffScalarSegregatedOp
{
public:
ICNSDiffScalarSegregatedOp(
PDEFields& fields,
const bool has_overset,
const bool mesh_mapping,
const std::string& prefix = "diffusion")
: m_pdefields(fields)
, m_density(fields.repo.get_field("density"))
, m_options(prefix, m_pdefields.field.name() + "_" + prefix)
, m_mesh_mapping(mesh_mapping)
{
amrex::LPInfo isolve = m_options.lpinfo();
amrex::LPInfo iapply;

iapply.setMaxCoarseningLevel(0);
const auto& mesh = m_pdefields.repo.mesh();

const auto& bclo = diffusion::get_diffuse_tensor_bc(
m_pdefields.field, amrex::Orientation::low);
const auto& bchi = diffusion::get_diffuse_tensor_bc(
m_pdefields.field, amrex::Orientation::high);

for (int i = 0; i < AMREX_SPACEDIM; ++i) {
if (!has_overset) {
m_solver_scalar[i].reset(new amrex::MLABecLaplacian(
mesh.Geom(0, mesh.finestLevel()),
mesh.boxArray(0, mesh.finestLevel()),
mesh.DistributionMap(0, mesh.finestLevel()), isolve));
m_applier_scalar[i].reset(new amrex::MLABecLaplacian(
mesh.Geom(0, mesh.finestLevel()),
mesh.boxArray(0, mesh.finestLevel()),
mesh.DistributionMap(0, mesh.finestLevel()), iapply));
} else {
auto imask =
fields.repo.get_int_field("mask_cell").vec_const_ptrs();
m_solver_scalar[i].reset(new amrex::MLABecLaplacian(
mesh.Geom(0, mesh.finestLevel()),
mesh.boxArray(0, mesh.finestLevel()),
mesh.DistributionMap(0, mesh.finestLevel()), imask,
isolve));
m_applier_scalar[i].reset(new amrex::MLABecLaplacian(
mesh.Geom(0, mesh.finestLevel()),
mesh.boxArray(0, mesh.finestLevel()),
mesh.DistributionMap(0, mesh.finestLevel()), imask,
iapply));
}

m_solver_scalar[i]->setMaxOrder(m_options.max_order);
m_applier_scalar[i]->setMaxOrder(m_options.max_order);

m_solver_scalar[i]->setDomainBC(bclo[i], bchi[i]);
m_applier_scalar[i]->setDomainBC(bclo[i], bchi[i]);
}
}

template <typename Scheme>
void compute_diff_term(const FieldState fstate)
{
auto tau_state = std::is_same<Scheme, fvm::Godunov>::value
? FieldState::New
: fstate;
const auto& repo = m_pdefields.repo;
const int nlevels = repo.num_active_levels();
const auto& geom = repo.mesh().Geom();

auto& divtau = m_pdefields.diff_term.state(tau_state);
const auto& density = m_density.state(fstate);
const auto& viscosity = m_pdefields.mueff;
Field const* mesh_detJ =
m_mesh_mapping ? &(repo.get_mesh_mapping_detJ(FieldLoc::CELL))
: nullptr;
std::unique_ptr<ScratchField> rho_times_detJ =
m_mesh_mapping ? repo.create_scratch_field(
1, m_density.num_grow()[0], FieldLoc::CELL)
: nullptr;

const amrex::Real alpha = 0.0;
const amrex::Real beta = -1.0;

for (int lev = 0; lev < nlevels; ++lev) {
if (m_mesh_mapping) {
(*rho_times_detJ)(lev).setVal(0.0);
amrex::MultiFab::AddProduct(
(*rho_times_detJ)(lev), density(lev), 0, (*mesh_detJ)(lev),
0, 0, 1, m_density.num_grow()[0]);
}
auto b = diffusion::average_velocity_eta_to_faces(
geom[lev], viscosity(lev));
if (m_mesh_mapping) {
diffusion::viscosity_to_uniform_space(b, repo, lev);
}

for (int i = 0; i < AMREX_SPACEDIM; ++i) {
m_applier_scalar[i]->setScalars(alpha, beta);

m_applier_scalar[i]->setLevelBC(
lev, &m_pdefields.field.subview(i)(lev));

// A coeffs
if (m_mesh_mapping) {
m_applier_scalar[i]->setACoeffs(
lev, (*rho_times_detJ)(lev));
} else {
m_applier_scalar[i]->setACoeffs(lev, density(lev));
}

// B coeffs
m_applier_scalar[i]->setBCoeffs(
lev, amrex::GetArrOfConstPtrs(b));

auto divtau_comp = divtau.subview(i);
auto vel_comp = m_pdefields.field.subview(i);

amrex::MLMG mlmg(*m_applier_scalar[i]);
mlmg.apply(divtau_comp.vec_ptrs(), vel_comp.vec_ptrs());
}
}

for (int lev = 0; lev < nlevels; ++lev) {
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for (amrex::MFIter mfi(divtau(lev), amrex::TilingIfNotGPU());
mfi.isValid(); ++mfi) {
const auto& bx = mfi.tilebox();
const auto& divtau_arr = divtau(lev).array(mfi);
const auto& rho_arr = density(lev).const_array(mfi);

amrex::ParallelFor(
bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
amrex::Real rhoinv = 1.0 / rho_arr(i, j, k);
divtau_arr(i, j, k, 0) *= rhoinv;
divtau_arr(i, j, k, 1) *= rhoinv;
divtau_arr(i, j, k, 2) *= rhoinv;
});
}
}
}

void linsys_solve(const amrex::Real dt)
{
const FieldState fstate = FieldState::New;
auto& repo = m_pdefields.repo;
auto& field = m_pdefields.field;
const auto& density = m_density.state(fstate);
const int nlevels = repo.num_active_levels();
const int ndim = field.num_comp();
auto rhs_ptr = repo.create_scratch_field("rhs", field.num_comp(), 0);
const auto& viscosity = m_pdefields.mueff;
const auto& geom = repo.mesh().Geom();

Field const* mesh_detJ =
m_mesh_mapping ? &(repo.get_mesh_mapping_detJ(FieldLoc::CELL))
: nullptr;
std::unique_ptr<ScratchField> rho_times_detJ =
m_mesh_mapping ? repo.create_scratch_field(
1, m_density.num_grow()[0], FieldLoc::CELL)
: nullptr;

const amrex::Real alpha = 1.0;
const amrex::Real beta = dt;

for (int lev = 0; lev < nlevels; ++lev) {

if (m_mesh_mapping) {
(*rho_times_detJ)(lev).setVal(0.0);
amrex::MultiFab::AddProduct(
(*rho_times_detJ)(lev), density(lev), 0, (*mesh_detJ)(lev),
0, 0, 1, m_density.num_grow()[0]);
}

auto b = diffusion::average_velocity_eta_to_faces(
geom[lev], viscosity(lev));
if (m_mesh_mapping) {
diffusion::viscosity_to_uniform_space(b, repo, lev);
}

for (int i = 0; i < AMREX_SPACEDIM; ++i) {
m_solver_scalar[i]->setScalars(alpha, beta);

m_solver_scalar[i]->setLevelBC(
lev, &m_pdefields.field.subview(i)(lev));

// A coeffs
if (m_mesh_mapping) {
m_solver_scalar[i]->setACoeffs(lev, (*rho_times_detJ)(lev));
} else {
m_solver_scalar[i]->setACoeffs(lev, density(lev));
}

// B coeffs
m_solver_scalar[i]->setBCoeffs(
lev, amrex::GetArrOfConstPtrs(b));
}
}

// Always multiply with rho since there is no diffusion term for density
for (int lev = 0; lev < nlevels; ++lev) {
auto& rhs = (*rhs_ptr)(lev);

#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for (amrex::MFIter mfi(rhs, amrex::TilingIfNotGPU()); mfi.isValid();
++mfi) {
const auto& bx = mfi.tilebox();
const auto& rhs_a = rhs.array(mfi);
const auto& fld = field(lev).const_array(mfi);
const auto& rho = density(lev).const_array(mfi);

amrex::ParallelFor(
bx, ndim,
[=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
rhs_a(i, j, k, n) = rho(i, j, k) * fld(i, j, k, n);
});
}
}

for (int i = 0; i < AMREX_SPACEDIM; ++i) {
auto vel_comp = m_pdefields.field.subview(i);
auto rhs_ptr_comp = rhs_ptr->subview(i);

amrex::MLMG mlmg(*m_solver_scalar[i]);
m_options(mlmg);

mlmg.solve(
vel_comp.vec_ptrs(), rhs_ptr_comp.vec_const_ptrs(),
m_options.rel_tol, m_options.abs_tol);

io::print_mlmg_info(
field.name() + std::to_string(i) + "_solve", mlmg);
}
}

protected:
PDEFields& m_pdefields;
Field& m_density;
MLMGOptions m_options;
bool m_mesh_mapping{false};

amrex::Array<std::unique_ptr<amrex::MLABecLaplacian>, AMREX_SPACEDIM>
m_solver_scalar;
amrex::Array<std::unique_ptr<amrex::MLABecLaplacian>, AMREX_SPACEDIM>
m_applier_scalar;
};

/** Specialization of diffusion operator for ICNS
* \ingroup icns
*/
Expand All @@ -293,25 +541,41 @@ struct DiffusionOp<ICNS, Scheme>
{
std::unique_ptr<ICNSDiffTensorOp> m_tensor_op;
std::unique_ptr<ICNSDiffScalarOp> m_scalar_op;
std::unique_ptr<ICNSDiffScalarSegregatedOp> m_scalar_segregated_op;

static_assert(
ICNS::ndim == AMREX_SPACEDIM,
"DiffusionOp invoked for scalar PDE type");

bool use_segregated_op = false;

DiffusionOp(
PDEFields& fields, const bool has_overset, const bool mesh_mapping)
{

bool use_tensor_op = true;

amrex::ParmParse pp(fields.field.name() + "_diffusion");
pp.query("use_tensor_operator", use_tensor_op);
pp.query("use_segregated_operator", use_segregated_op);

if (use_tensor_op && use_segregated_op) {
amrex::Abort(
"Tensor and segregated operators should not be enabled "
"simultaneously.");
}

if (use_tensor_op) {
m_tensor_op = std::make_unique<ICNSDiffTensorOp>(
fields, has_overset, mesh_mapping);
} else {
m_scalar_op = std::make_unique<ICNSDiffScalarOp>(
fields, has_overset, mesh_mapping);
if (use_segregated_op) {
m_scalar_segregated_op =
std::make_unique<ICNSDiffScalarSegregatedOp>(
fields, has_overset, mesh_mapping);
} else {
m_scalar_op = std::make_unique<ICNSDiffScalarOp>(
fields, has_overset, mesh_mapping);
}
}
}

Expand All @@ -320,7 +584,11 @@ struct DiffusionOp<ICNS, Scheme>
if (m_tensor_op) {
m_tensor_op->compute_diff_term<Scheme>(fstate);
} else {
m_scalar_op->compute_diff_term<Scheme>(fstate);
if (use_segregated_op) {
m_scalar_segregated_op->compute_diff_term<Scheme>(fstate);
} else {
m_scalar_op->compute_diff_term<Scheme>(fstate);
}
}
}

Expand All @@ -329,7 +597,11 @@ struct DiffusionOp<ICNS, Scheme>
if (m_tensor_op) {
m_tensor_op->linsys_solve(dt);
} else {
m_scalar_op->linsys_solve(dt);
if (use_segregated_op) {
m_scalar_segregated_op->linsys_solve(dt);
} else {
m_scalar_op->linsys_solve(dt);
}
}
}
};
Expand Down
2 changes: 2 additions & 0 deletions test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,7 @@ add_test_re(abl_godunov_nolim)
add_test_re(abl_godunov_plm)
add_test_re(abl_godunov_static_refinement)
add_test_re(abl_godunov_scalar_velocity_solve)
add_test_re(abl_godunov_segregated_velocity_solve)
add_test_re(abl_ksgsm84_godunov)
add_test_re(abl_mol_cn)
add_test_re(abl_mol_explicit)
Expand Down Expand Up @@ -224,6 +225,7 @@ if(AMR_WIND_ENABLE_HYPRE)
add_test_re(channel_mol_mesh_map_x)
add_test_re(channel_mol_mesh_map_y)
add_test_re(channel_mol_mesh_map_z)
add_test_re(channel_mol_mesh_map_x_seg_vel_solve)
endif()

if(AMR_WIND_ENABLE_HDF5)
Expand Down
Loading

0 comments on commit 7b68af2

Please sign in to comment.