8#ifndef AMREX_HYPRE_SOLVER_H_
9#define AMREX_HYPRE_SOLVER_H_
17#include "_hypre_utilities.h"
63 template <
class Marker,
class Filler>
72 std::string a_options_namespace =
"hypre");
74 template <class MF, std::enable_if_t<IsFabArray<MF>::value &&
75 std::is_same_v<
typename MF::value_type,
76 HYPRE_Real>,
int> = 0>
89 HYPRE_Real rel_tol, HYPRE_Real abs_tol,
int max_iter);
103 template <
class Marker>
105 std::enable_if_t<IsCallable<Marker,int,int,int,int,int>::value>
107 std::enable_if_t<IsCallableR<bool,Marker,int,int,int,int,int>::value>
111 template <
typename AI>
114 template <
class Filler,
117 HYPRE_Int&, HYPRE_Int*,
118 HYPRE_Real*>::value,
int> FOO = 0>
121 template <class MF, std::enable_if_t<IsFabArray<MF>::value &&
122 std::is_same_v<
typename MF::value_type,
123 HYPRE_Real>,
int> = 0>
127 template <class MF, std::enable_if_t<IsFabArray<MF>::value &&
128 std::is_same_v<
typename MF::value_type,
129 HYPRE_Real>,
int> = 0>
163 HYPRE_IJMatrix
m_A =
nullptr;
164 HYPRE_IJVector
m_b =
nullptr;
165 HYPRE_IJVector
m_x =
nullptr;
169template <
class Marker,
class Filler>
178 std::string a_options_namespace)
179 : m_nvars (
int(a_index_type.size())),
180 m_index_type (a_index_type),
184 m_verbose (a_verbose),
185 m_options_namespace(std::move(a_options_namespace))
199 for (
int ivar = 0; ivar <
m_nvars; ++ivar) {
205 nrows_max +=
m_grids[ivar].numPts();
210 "Need to configure Hypre with --enable-bigint");
213 for (
int ivar = 0; ivar <
m_nvars; ++ivar) {
235 if (nrows_allprocs.
size() > 1) {
236 MPI_Allgather(&
m_nrows_proc,
sizeof(HYPRE_Int), MPI_CHAR,
237 nrows_allprocs.data(),
sizeof(HYPRE_Int), MPI_CHAR,
m_comm);
244 HYPRE_Int proc_begin = 0;
245 for (
int i = 0; i < myproc; ++i) {
246 proc_begin += nrows_allprocs[i];
249 HYPRE_Int proc_end = proc_begin;
251 for (
int ivar = 0; ivar <
m_nvars; ++ivar) {
263 using AtomicInt = std::conditional_t<
sizeof(HYPRE_Int) == 4,
264 HYPRE_Int,
unsigned long long>;
265 fill_global_id<AtomicInt>();
268 HYPRE_Int ilower = proc_begin;
269 HYPRE_Int iupper = proc_end-1;
282template <
class Marker>
284 std::enable_if_t<IsCallable<Marker,int,int,int,int,int>::value>
286 std::enable_if_t<IsCallableR<bool,Marker,int,int,int,int,int>::value>
295 int boxno = mfi.LocalIndex();
297 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
301 m_cell_offset[mfi].resize(npts_tot);
303 int* p_cell_offset = m_cell_offset[mfi].data();
304 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
306 auto const& lid = m_local_id[ivar].array(mfi);
307 auto const& owner = m_owner_mask[ivar]->const_array(mfi);
309 const auto npts =
static_cast<int>(bx.
numPts());
310 int npts_box = amrex::Scan::PrefixSum<int>(npts,
314 int id = (owner ( cell.
x,cell.
y,cell.
z ) &&
315 marker(boxno,cell.
x,cell.
y,cell.
z,ivar)) ? 1 : 0;
316 lid(cell.
x,cell.
y,cell.
z) = id;
322 if (lid(cell.
x,cell.
y,cell.
z)) {
323 lid(cell.
x,cell.
y,cell.
z) = ps;
324 p_cell_offset[ps] =
offset;
326 lid(cell.
x,cell.
y,cell.
z) = std::numeric_limits<int>::lowest();
330 m_nrows_grid[ivar][mfi] = npts_box;
331 npts_tot += npts_box;
332 p_cell_offset += npts_box;
334 m_cell_offset[mfi].resize(npts_tot);
343 int boxno = mfi.LocalIndex();
344 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
346 auto const& lid = m_local_id[ivar].array(mfi);
347 auto const& owner = m_owner_mask[ivar]->const_array(mfi);
351 for (
int k = lo.z; k <= hi.z; ++k) {
352 for (
int j = lo.y; j <= hi.y; ++j) {
353 for (
int i = lo.x; i <= hi.x; ++i) {
354 if (owner(i,j,k) && marker(boxno,i,j,k,ivar)) {
357 lid(i,j,k) = std::numeric_limits<int>::lowest();
360 m_nrows_grid[ivar][mfi] = id;
368 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
369 nrows += m_nrows_grid[ivar][mfi];
371 m_nrows[mfi] = nrows;
372 m_nrows_proc += nrows;
377template <
typename AI>
386 if constexpr (std::is_same_v<HYPRE_Int,AI>) {
387 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
388 p_global_id.push_back(&(m_global_id[ivar]));
391 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
392 global_id_raii.emplace_back(m_global_id[ivar].
boxArray(),
393 m_global_id[ivar].DistributionMap(),
394 1, m_global_id[ivar].nGrowVect());
395 p_global_id.push_back(&(global_id_raii[ivar]));
400#pragma omp parallel if (Gpu::notInLaunchRegion())
403 auto& rows_vec = m_global_id_vec[mfi];
404 rows_vec.resize(m_nrows[mfi]);
407 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
408 HYPRE_Int
const os = m_id_offset[ivar][mfi];
409 Box bx = mfi.validbox();
411 Array4<AI> const& gid = p_global_id[ivar]->array(mfi);
412 auto const& lid = m_local_id[ivar].const_array(mfi);
413 HYPRE_Int* rows = rows_vec.data() + nrows;
414 nrows += m_nrows_grid[ivar][mfi];
417 if (lid.contains(i,j,k) && lid(i,j,k) >= 0) {
418 const auto id = lid(i,j,k) + os;
419 rows[lid(i,j,k)] = id;
420 gid(i,j,k) =
static_cast<AI
>(id);
422 gid(i,j,k) =
static_cast<AI
>
423 (std::numeric_limits<HYPRE_Int>::max());
429 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
431 m_geom.periodicity());
432 p_global_id[ivar]->FillBoundary(m_geom.periodicity());
434 if constexpr (!std::is_same<HYPRE_Int, AI>()) {
435 auto const& dst = m_global_id[ivar].arrays();
436 auto const& src = p_global_id[ivar]->const_arrays();
440 dst[
b](i,j,k) =
static_cast<HYPRE_Int
>(src[
b](i,j,k));
454 auto* p_cols_tmp = cols_tmp.
data();
455 auto* p_mat_tmp = mat_tmp.
data();
456 auto const* p_cols = cols.
data();
457 auto const* p_mat = mat.
data();
458 const auto N = Long(cols.
size());
459 Scan::PrefixSum<T>(N,
462 return static_cast<T
>(p_cols[i] >= 0);
466 if (p_cols[i] >= 0) {
467 p_cols_tmp[s] = p_cols[i];
468 p_mat_tmp[s] = p_mat[i];
472 std::swap(cols_tmp, cols);
473 std::swap(mat_tmp, mat);
479template <
class Filler,
481 Array4<HYPRE_Int const>
const*,
482 HYPRE_Int&, HYPRE_Int*,
483 HYPRE_Real*>::value,
int> FOO>
495 for (
MFIter mfi(m_local_id[0],mfitinfo); mfi.
isValid(); ++mfi)
497 int boxno = mfi.LocalIndex();
498 const HYPRE_Int nrows = m_nrows[mfi];
503 HYPRE_Int* ncols = ncols_vec.
data();
506 cols_vec.
resize(Long(nrows)*MSS, -1);
507 HYPRE_Int* cols = cols_vec.
data();
510 mat_vec.
resize(Long(nrows)*MSS);
511 HYPRE_Real* mat = mat_vec.
data();
514 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
515 gid_v[ivar] = m_global_id[ivar].const_array(mfi);
520 (gid_v.data(), gid_v.
size());
521 auto const* pgid = gid_buf.
data();
522 auto const* p_cell_offset = m_cell_offset[mfi].data();
524 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
525 const HYPRE_Int nrows_var = m_nrows_grid[ivar][mfi];
528 ntot += Reduce::Sum<Long>(nrows_var,
532 filler(boxno, cell.
x, cell.
y, cell.
z, ivar, pgid,
537 p_cell_offset += nrows_var;
539 cols += Long(nrows_var)*MSS;
540 mat += Long(nrows_var)*MSS;
545 if (ntot >= Long(std::numeric_limits<int>::max())) {
546 detail::pack_matrix_gpu<Long>(cols_tmp, mat_tmp, cols_vec, mat_vec);
548 detail::pack_matrix_gpu<int>(cols_tmp, mat_tmp, cols_vec, mat_vec);
551 auto* pgid = gid_v.data();
552 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
553 if (m_nrows_grid[ivar][mfi] > 0) {
554 auto const& lid = m_local_id[ivar].const_array(mfi);
556 [=,&ncols,&cols,&mat] (
int i,
int j,
int k)
558 if (lid(i,j,k) >= 0) {
559 filler(boxno, i, j, k, ivar, pgid, *ncols, cols, mat);
569 const auto& rows_vec = m_global_id_vec[mfi];
570 HYPRE_Int
const* rows = rows_vec.data();
573 HYPRE_IJMatrixSetValues(m_A, nrows, ncols_vec.
data(), rows,
575 Gpu::hypreSynchronize();
578 HYPRE_IJMatrixAssemble(m_A);
582template <class MF, std::enable_if_t<IsFabArray<MF>::value &&
583 std::is_same_v<
typename MF::value_type,
584 HYPRE_Real>,
int> FOO>
588 HYPRE_Real rel_tol, HYPRE_Real abs_tol,
int max_iter)
594 HYPRE_IJVectorInitialize(m_b);
595 HYPRE_IJVectorInitialize(m_x);
597 load_vectors(a_soln, a_rhs);
599 HYPRE_IJVectorAssemble(m_x);
600 HYPRE_IJVectorAssemble(m_b);
602 m_hypre_ij->solve(rel_tol, abs_tol, max_iter);
604 get_solution(a_soln);
608template <class MF, std::enable_if_t<IsFabArray<MF>::value &&
609 std::is_same_v<
typename MF::value_type,
610 HYPRE_Real>,
int> FOO>
624 const HYPRE_Int nrows = m_nrows[mfi];
631 auto* xp = xvec.
data();
632 auto* bp = bvec.
data();
634 HYPRE_Int
const* rows = m_global_id_vec[mfi].data();
637 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
638 if (m_nrows_grid[ivar][mfi] > 0) {
639 auto const& xfab = a_soln[ivar]->const_array(mfi);
640 auto const& bfab = a_rhs [ivar]->const_array(mfi);
641 auto const& lid = m_local_id[ivar].const_array(mfi);
647 if (lid(i,j,k) >= 0) {
648 x[lid(i,j,k)] = xfab(i,j,k);
649 b[lid(i,j,k)] = bfab(i,j,k);
652 offset += m_nrows_grid[ivar][mfi];
656 Gpu::streamSynchronize();
657 HYPRE_IJVectorSetValues(m_x, nrows, rows, xp);
658 HYPRE_IJVectorSetValues(m_b, nrows, rows, bp);
659 Gpu::hypreSynchronize();
665template <class MF, std::enable_if_t<IsFabArray<MF>::value &&
666 std::is_same_v<
typename MF::value_type,
667 HYPRE_Real>,
int> FOO>
679 const HYPRE_Int nrows = m_nrows[mfi];
684 auto* xp = xvec.
data();
686 HYPRE_Int
const* rows = m_global_id_vec[mfi].data();
688 HYPRE_IJVectorGetValues(m_x, nrows, rows, xp);
689 Gpu::hypreSynchronize();
692 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
693 if (m_nrows_grid[ivar][mfi] > 0) {
694 auto const& xfab = a_soln[ivar]->array(mfi);
695 auto const& lid = m_local_id[ivar].const_array(mfi);
700 if (lid(i,j,k) >= 0) {
701 xfab(i,j,k) =
x[lid(i,j,k)];
704 offset += m_nrows_grid[ivar][mfi];
707 Gpu::streamSynchronize();
711 for (
int ivar = 0; ivar < m_nvars; ++ivar) {
713 m_geom.periodicity());
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:49
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
int MPI_Comm
Definition AMReX_ccse-mpi.H:47
static constexpr int MPI_COMM_NULL
Definition AMReX_ccse-mpi.H:55
A collection of Boxes stored in an Array.
Definition AMReX_BoxArray.H:550
AMREX_GPU_HOST_DEVICE IntVectND< dim > atOffset(Long offset) const noexcept
Given the offset, compute IntVectND<dim>
Definition AMReX_Box.H:1009
AMREX_GPU_HOST_DEVICE BoxND & convert(IndexTypeND< dim > typ) noexcept
Convert the BoxND from the current type into the argument type. This may change the BoxND coordinates...
Definition AMReX_Box.H:912
AMREX_GPU_HOST_DEVICE BoxND & grow(int i) noexcept
Definition AMReX_Box.H:627
AMREX_GPU_HOST_DEVICE Long numPts() const noexcept
Returns the number of points contained in the BoxND.
Definition AMReX_Box.H:346
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:41
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:73
Periodicity periodicity() const noexcept
Definition AMReX_Geometry.H:355
Definition AMReX_GpuBuffer.H:17
T const * data() const noexcept
Definition AMReX_GpuBuffer.H:65
Solve Ax = b using HYPRE's generic IJ matrix format where A is a sparse matrix specified using the co...
Definition AMReX_HypreSolver.H:34
HypreSolver(Vector< IndexType > const &a_index_type, IntVect const &a_nghost, Geometry const &a_geom, BoxArray const &a_grids, DistributionMapping const &a_dmap, Marker &&a_marker, Filler &&a_filler, int a_verbose=0, std::string a_options_namespace="hypre")
Definition AMReX_HypreSolver.H:170
LayoutData< HYPRE_Int > m_nrows
Definition AMReX_HypreSolver.H:157
IntVect m_nghost
Definition AMReX_HypreSolver.H:136
Vector< LayoutData< HYPRE_Int > > m_nrows_grid
Definition AMReX_HypreSolver.H:155
int getNumIters() const
Definition AMReX_HypreSolver.H:91
Geometry m_geom
Definition AMReX_HypreSolver.H:137
LayoutData< Gpu::DeviceVector< HYPRE_Int > > m_global_id_vec
Definition AMReX_HypreSolver.H:149
std::enable_if_t< IsCallable< Marker, int, int, int, int, int >::value > fill_local_id(Marker const &marker)
Definition AMReX_HypreSolver.H:288
std::unique_ptr< HypreIJIface > m_hypre_ij
Definition AMReX_HypreSolver.H:160
int m_nvars
Definition AMReX_HypreSolver.H:134
Vector< LayoutData< HYPRE_Int > > m_id_offset
Definition AMReX_HypreSolver.H:156
MPI_Comm m_comm
Definition AMReX_HypreSolver.H:144
void get_solution(Vector< MF * > const &a_soln)
Definition AMReX_HypreSolver.H:669
HYPRE_IJVector getx() const
Definition AMReX_HypreSolver.H:99
void load_vectors(Vector< MF * > const &a_soln, Vector< MF const * > const &a_rhs)
Definition AMReX_HypreSolver.H:612
HYPRE_IJMatrix getA() const
Definition AMReX_HypreSolver.H:97
std::string m_options_namespace
Definition AMReX_HypreSolver.H:142
void fill_matrix(Filler const &filler)
Definition AMReX_HypreSolver.H:485
HYPRE_IJMatrix m_A
Definition AMReX_HypreSolver.H:163
void fill_global_id()
Definition AMReX_HypreSolver.H:379
void solve(Vector< MF * > const &a_soln, Vector< MF const * > const &a_rhs, HYPRE_Real rel_tol, HYPRE_Real abs_tol, int max_iter)
Definition AMReX_HypreSolver.H:586
HYPRE_Int m_nrows_proc
Definition AMReX_HypreSolver.H:158
LayoutData< Gpu::DeviceVector< int > > m_cell_offset
Definition AMReX_HypreSolver.H:152
Vector< FabArray< BaseFab< HYPRE_Int > > > m_global_id
Definition AMReX_HypreSolver.H:148
HYPRE_Real getFinalResidualNorm() const
Definition AMReX_HypreSolver.H:93
Vector< std::unique_ptr< iMultiFab > > m_owner_mask
Definition AMReX_HypreSolver.H:146
HYPRE_IJVector m_b
Definition AMReX_HypreSolver.H:164
Vector< IndexType > m_index_type
Definition AMReX_HypreSolver.H:135
Vector< BoxArray > m_grids
Definition AMReX_HypreSolver.H:138
Vector< iMultiFab > m_local_id
Definition AMReX_HypreSolver.H:147
int m_verbose
Definition AMReX_HypreSolver.H:141
DistributionMapping m_dmap
Definition AMReX_HypreSolver.H:139
HYPRE_IJVector m_x
Definition AMReX_HypreSolver.H:165
HYPRE_IJVector getb() const
Definition AMReX_HypreSolver.H:98
a one-thingy-per-box distributed object
Definition AMReX_LayoutData.H:13
void define(const BoxArray &a_grids, const DistributionMapping &a_dm)
Definition AMReX_LayoutData.H:25
Definition AMReX_MFIter.H:57
bool isValid() const noexcept
Is the iterator valid i.e. is it associated with a FAB?
Definition AMReX_MFIter.H:141
Definition AMReX_PODVector.H:262
size_type size() const noexcept
Definition AMReX_PODVector.H:591
void clear() noexcept
Definition AMReX_PODVector.H:589
void resize(size_type a_new_size)
Definition AMReX_PODVector.H:641
T * data() noexcept
Definition AMReX_PODVector.H:609
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:27
Long size() const noexcept
Definition AMReX_Vector.H:50
void streamSynchronize() noexcept
Definition AMReX_GpuDevice.H:237
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition AMReX_ParallelContext.H:70
int MyProcSub() noexcept
my sub-rank in current frame
Definition AMReX_ParallelContext.H:76
int NProcsSub() noexcept
number of ranks in current frame
Definition AMReX_ParallelContext.H:74
static constexpr struct amrex::Scan::Type::Exclusive exclusive
static constexpr RetSum noRetSum
Definition AMReX_Scan.H:30
void pack_matrix_gpu(Gpu::DeviceVector< HYPRE_Int > &cols_tmp, Gpu::DeviceVector< HYPRE_Real > mat_tmp, Gpu::DeviceVector< HYPRE_Int > &cols, Gpu::DeviceVector< HYPRE_Real > &mat)
Definition AMReX_HypreSolver.H:449
Definition AMReX_Amr.cpp:49
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
Definition AMReX_CTOParallelForImpl.H:191
IntVect nGrowVect(FabArrayBase const &fa)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Returns a BoxND with different type.
Definition AMReX_Box.H:1435
void OverrideSync(FabArray< FAB > &fa, FabArray< IFAB > const &msk, const Periodicity &period)
Definition AMReX_FabArrayUtility.H:1323
std::unique_ptr< iMultiFab > OwnerMask(FabArrayBase const &mf, const Periodicity &period, const IntVect &ngrow)
Definition AMReX_iMultiFab.cpp:637
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:308
const int[]
Definition AMReX_BLProfiler.cpp:1664
BoxArray const & boxArray(FabArrayBase const &fa)
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void Loop(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:125
Definition AMReX_FabArrayCommI.H:896
Definition AMReX_Array4.H:61
Definition AMReX_Dim3.H:12
int x
Definition AMReX_Dim3.H:12
int z
Definition AMReX_Dim3.H:12
int y
Definition AMReX_Dim3.H:12
Test if a given type T is callable with arguments of type Args...
Definition AMReX_TypeTraits.H:201
Definition AMReX_MFIter.H:20
MFItInfo & DisableDeviceSync() noexcept
Definition AMReX_MFIter.H:38
MFItInfo & UseDefaultStream() noexcept
Definition AMReX_MFIter.H:50