1 #ifndef AMREX_EB2_GEOMETRYSHOP_H_
2 #define AMREX_EB2_GEOMETRYSHOP_H_
3 #include <AMReX_Config.H>
12 #include <type_traits>
18 template <class F, std::enable_if_t<IsGPUable<F>::value>* FOO =
nullptr>
26 template <class F, std::enable_if_t<!IsGPUable<F>::value>* BAR =
nullptr>
44 int rangedir, F
const&
f) noexcept
46 #ifdef AMREX_USE_FLOAT
47 const Real tol = 1.e-4_rt;
48 const Real EPS = 1.0e-6_rt;
50 const Real tol = 1.e-12;
51 const Real EPS = 3.0e-15;
53 const int MAXITER = 100;
60 Real fa =
IF_f(
f, aPt);
61 Real fb =
IF_f(
f, bPt);
62 Real c = bPt[rangedir];
67 amrex::Error(
"BrentRootFinder. Root must be bracketed, but instead the supplied end points have the same sign.");
69 }
else if (fa == 0.0_rt) {
71 }
else if (fb == 0.0_rt) {
75 Real d = 0.0_rt, e = 0.0_rt;
77 for (i = 0; i < MAXITER; ++i)
84 d = bPt[rangedir] - aPt[rangedir];
90 aPt[rangedir] = bPt[rangedir];
99 Real tol1 = 2.0_rt * EPS *
std::abs(bPt[rangedir]) + 0.5_rt * tol;
100 Real xm = 0.5_rt * (c - bPt[rangedir]);
102 if (
std::abs(xm) <= tol1 || fb == 0.0_rt)
111 if (aPt[rangedir] == c)
120 p = s * (2.0_rt * xm * q * (q-
r) - (bPt[rangedir]-aPt[rangedir]) * (
r-1.0_rt));
121 q = (q-1.0_rt) * (
r-1.0_rt) * (s-1.0_rt);
125 if (p > 0) { q = -q; }
150 aPt[rangedir] = bPt[rangedir];
156 bPt[rangedir] = bPt[rangedir] + d;
160 if (xm < 0) { bPt[rangedir] = bPt[rangedir] - tol1; }
161 else { bPt[rangedir] = bPt[rangedir] + tol1; }
169 amrex::Error(
"BrentRootFinder: exceeding maximum iterations.");
172 return bPt[rangedir];
175 template <
class F,
class R =
int>
203 const Real* problo = geom.ProbLo();
204 const Real* dx = geom.CellSize();
205 const auto& len3 = bx.length3d();
206 const int* blo = bx.loVect();
207 int nbody = 0, nzero = 0, nfluid = 0;
208 for (
int k = 0; k < len3[2]; ++k) {
209 for (
int j = 0; j < len3[1]; ++j) {
210 for (
int i = 0; i < len3[0]; ++i) {
212 problo[1]+(j+blo[1])*dx[1],
213 problo[2]+(k+blo[2])*dx[2])};
217 }
else if (v > 0.0_rt) {
222 if (nbody > 0 && nfluid > 0) {
return mixedcells; }
230 }
else if (nfluid == 0) {
237 template <class U=F, std::enable_if_t<IsGPUable<U>::value>* FOO =
nullptr >
242 const auto& problo = geom.ProbLoArray();
243 const auto& dx = geom.CellSizeArray();
247 using ReduceTuple =
typename decltype(reduce_data)::Type;
248 reduce_op.
eval(bx, reduce_data,
255 int nbody = 0, nzero = 0, nfluid = 0;
258 }
else if (v > 0.0_rt) {
263 return {nbody,nzero,nfluid};
265 ReduceTuple hv = reduce_data.
value(reduce_op);
266 int nbody = amrex::get<0>(hv);
268 int nfluid = amrex::get<2>(hv);
272 }
else if (nfluid == 0) {
284 template <class U=F, std::enable_if_t<!IsGPUable<U>::value>* BAR =
nullptr >
290 template <class U=F, std::enable_if_t<IsGPUable<U>::value>* FOO =
nullptr >
291 static constexpr
bool isGPUable () noexcept {
return true; }
293 template <class U=F, std::enable_if_t<!IsGPUable<U>::value>* BAR =
nullptr >
294 static constexpr
bool isGPUable () noexcept {
return false; }
296 template <class U=F, std::enable_if_t<IsGPUable<U>::value>* FOO =
nullptr >
298 Box const& bounding_box)
const noexcept
300 const auto& problo = geom.ProbLoArray();
301 const auto& dx = geom.CellSizeArray();
302 const Box& bx = levelset.box();
303 const auto& a = levelset.array();
315 template <class U=F, std::enable_if_t<!IsGPUable<U>::value>* BAR =
nullptr >
317 Box const& bounding_box)
const noexcept
320 if (!levelset.arena()->isHostAccessible()) {
321 const Box& bx = levelset.box();
325 levelset.nBytes(bx, levelset.nComp()));
336 Box const& bounding_box)
const noexcept
338 const auto& problo = geom.ProbLoArray();
339 const auto& dx = geom.CellSizeArray();
340 const Box& bx = levelset.box();
344 const auto& a = levelset.array();
354 template <class U=F, std::enable_if_t<IsGPUable<U>::value>* FOO =
nullptr >
358 Box const& bounding_box)
const noexcept
360 auto const& dx = geom.CellSizeArray();
361 auto const& problo = geom.ProbLoArray();
365 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
384 inter(i,j,k) = std::numeric_limits<Real>::quiet_NaN();
390 template <class U=F, std::enable_if_t<!IsGPUable<U>::value>* BAR =
nullptr >
394 Box const& bounding_box)
const noexcept
396 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
407 h_type_fab.
nBytes(bx, 1));
409 const auto& h_inter = h_inter_fab.
array();
413 h_inter_fab.
nBytes(bx, 1));
428 Box const& bounding_box,
429 const int idim)
const noexcept
431 auto const& dx = geom.CellSizeArray();
432 auto const& problo = geom.ProbLoArray();
451 inter(i,j,k) = std::numeric_limits<Real>::quiet_NaN();
462 auto const& dx = geom.CellSizeArray();
463 auto const& problo = geom.ProbLoArray();
464 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
473 if (lst(i,j,k) == Real(0.0) ||
474 (lst(i,j,k) > Real(0.0) && is_nan))
479 inter(i,j,k) = problo[0] + i*dx[0];
481 else if (lst(i+1,j,k) == Real(0.0) ||
482 (lst(i+1,j,k) > Real(0.0) && is_nan))
484 inter(i,j,k) = problo[0] + (i+1)*dx[0];
486 }
else if (idim == 1) {
487 if (lst(i,j,k) == Real(0.0) ||
488 (lst(i,j,k) > Real(0.0) && is_nan))
490 inter(i,j,k) = problo[1] + j*dx[1];
492 else if (lst(i,j+1,k) == Real(0.0) ||
493 (lst(i,j+1,k) > Real(0.0) && is_nan))
495 inter(i,j,k) = problo[1] + (j+1)*dx[1];
498 if (lst(i,j,k) == Real(0.0) ||
499 (lst(i,j,k) > Real(0.0) && is_nan))
501 inter(i,j,k) = problo[2] + k*dx[2];
503 else if (lst(i,j,k+1) == Real(0.0) ||
504 (lst(i,j,k+1) > Real(0.0) && is_nan))
506 inter(i,j,k) = problo[2] + (k+1)*dx[2];
528 template <
class F,
class R>
529 GeometryShop<std::decay_t<F>, std::decay_t<R>>
534 (std::forward<F>(
f), std::forward<R>(
r));
#define AMREX_HOST_DEVICE_FOR_3D_FLAG(where_to_run, box, i, j, k, block)
Definition: AMReX_GpuLaunch.nolint.H:116
#define AMREX_HOST_DEVICE_PARALLEL_FOR_3D(...)
Definition: AMReX_GpuLaunch.nolint.H:54
#define AMREX_IF_ON_DEVICE(CODE)
Definition: AMReX_GpuQualifiers.H:56
#define AMREX_IF_ON_HOST(CODE)
Definition: AMReX_GpuQualifiers.H:58
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
#define AMREX_D_DECL(a, b, c)
Definition: AMReX_SPACE.H:104
AMREX_FORCE_INLINE Array4< T const > const_array() const noexcept
Definition: AMReX_BaseFab.H:415
std::size_t nBytes() const noexcept
Returns how many bytes used.
Definition: AMReX_BaseFab.H:266
T * dataPtr(int n=0) noexcept
Returns a pointer to an object of type T that is the value of the Nth component associated with the c...
Definition: AMReX_BaseFab.H:352
AMREX_FORCE_INLINE Array4< T const > array() const noexcept
Definition: AMReX_BaseFab.H:379
Definition: AMReX_EB2_GeometryShop.H:177
F && GetImpFunc() &&
Definition: AMReX_EB2_GeometryShop.H:199
void fillFab(BaseFab< Real > &levelset, const Geometry &geom, RunOn, Box const &bounding_box) const noexcept
Definition: AMReX_EB2_GeometryShop.H:316
static constexpr bool isGPUable() noexcept
Definition: AMReX_EB2_GeometryShop.H:291
GeometryShop(F f, R r)
Definition: AMReX_EB2_GeometryShop.H:194
static constexpr int in_body
Definition: AMReX_EB2_GeometryShop.H:182
void getIntercept(Array< Array4< Real >, AMREX_SPACEDIM > const &inter_arr, Array< Array4< Type_t const >, AMREX_SPACEDIM > const &type_arr, Array4< Real const > const &, Geometry const &geom, RunOn, Box const &bounding_box) const noexcept
Definition: AMReX_EB2_GeometryShop.H:391
static constexpr int allcovered
Definition: AMReX_EB2_GeometryShop.H:186
int getBoxType(const Box &bx, const Geometry &geom, RunOn) const noexcept
Definition: AMReX_EB2_GeometryShop.H:285
void updateIntercept(Array< Array4< Real >, AMREX_SPACEDIM > const &inter_arr, Array< Array4< Type_t const >, AMREX_SPACEDIM > const &type_arr, Array4< Real const > const &lst, Geometry const &geom) const noexcept
Definition: AMReX_EB2_GeometryShop.H:458
GeometryShop(F f)
Definition: AMReX_EB2_GeometryShop.H:190
static constexpr int in_fluid
Definition: AMReX_EB2_GeometryShop.H:180
R m_resource
Definition: AMReX_EB2_GeometryShop.H:517
void fillFab_Cpu(BaseFab< Real > &levelset, const Geometry &geom, Box const &bounding_box) const noexcept
Definition: AMReX_EB2_GeometryShop.H:335
F const & GetImpFunc() const &
Definition: AMReX_EB2_GeometryShop.H:198
F FunctionType
Definition: AMReX_EB2_GeometryShop.H:188
void fillFab(BaseFab< Real > &levelset, const Geometry &geom, RunOn run_on, Box const &bounding_box) const noexcept
Definition: AMReX_EB2_GeometryShop.H:297
static constexpr int mixedcells
Definition: AMReX_EB2_GeometryShop.H:185
static constexpr int allregular
Definition: AMReX_EB2_GeometryShop.H:184
void getIntercept_Cpu(Array4< Real > const &inter, Array4< Type_t const > const &type, Geometry const &geom, Box const &bounding_box, const int idim) const noexcept
Definition: AMReX_EB2_GeometryShop.H:425
static constexpr int on_boundary
Definition: AMReX_EB2_GeometryShop.H:181
void getIntercept(Array< Array4< Real >, AMREX_SPACEDIM > const &inter_arr, Array< Array4< Type_t const >, AMREX_SPACEDIM > const &type_arr, Array4< Real const > const &, Geometry const &geom, RunOn run_on, Box const &bounding_box) const noexcept
Definition: AMReX_EB2_GeometryShop.H:355
int getBoxType_Cpu(const Box &bx, Geometry const &geom) const noexcept
Definition: AMReX_EB2_GeometryShop.H:201
F m_f
Definition: AMReX_EB2_GeometryShop.H:516
int getBoxType(const Box &bx, const Geometry &geom, RunOn run_on) const noexcept
Definition: AMReX_EB2_GeometryShop.H:238
Rectangular problem domain geometry.
Definition: AMReX_Geometry.H:73
Definition: AMReX_Reduce.H:249
Type value()
Definition: AMReX_Reduce.H:281
Definition: AMReX_Reduce.H:364
std::enable_if_t< IsFabArray< MF >::value > eval(MF const &mf, IntVect const &nghost, D &reduce_data, F &&f)
Definition: AMReX_Reduce.H:441
static constexpr Type_t irregular
Definition: AMReX_EB2_Graph.H:40
Definition: AMReX_FabArrayBase.H:32
AMREX_GPU_HOST_DEVICE Real IF_f(F const &f, GpuArray< Real, AMREX_SPACEDIM > const &p) noexcept
Definition: AMReX_EB2_GeometryShop.H:21
GeometryShop< std::decay_t< F > > makeShop(F &&f)
Definition: AMReX_EB2_GeometryShop.H:523
AMREX_GPU_HOST_DEVICE Real BrentRootFinder(GpuArray< Real, AMREX_SPACEDIM > const &lo, GpuArray< Real, AMREX_SPACEDIM > const &hi, int rangedir, F const &f) noexcept
Definition: AMReX_EB2_GeometryShop.H:42
void streamSynchronize() noexcept
Definition: AMReX_GpuDevice.H:237
void dtoh_memcpy_async(void *p_h, const void *p_d, const std::size_t sz) noexcept
Definition: AMReX_GpuDevice.H:265
bool inLaunchRegion() noexcept
Definition: AMReX_GpuControl.H:86
void htod_memcpy_async(void *p_d, const void *p_h, const std::size_t sz) noexcept
Definition: AMReX_GpuDevice.H:251
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isnan(T m) noexcept
Definition: AMReX_GpuUtility.H:150
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
AMREX_ATTRIBUTE_FLATTEN_FOR void LoopOnCpu(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:355
RunOn
Definition: AMReX_GpuControl.H:69
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition: AMReX_GpuComplex.H:356
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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition: AMReX.H:111
Array< Real, AMREX_SPACEDIM > RealArray
Definition: AMReX_Array.H:26
void Error(const std::string &msg)
Print out message to cerr and exit via amrex::Abort().
Definition: AMReX.cpp:219
Arena * The_Pinned_Arena()
Definition: AMReX_Arena.cpp:649
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & Clamp(const T &v, const T &lo, const T &hi)
Definition: AMReX_Algorithm.H:84
Arena * The_Arena()
Definition: AMReX_Arena.cpp:609
std::array< T, N > Array
Definition: AMReX_Array.H:24
Definition: AMReX_Array4.H:61
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T * dataPtr() const noexcept
Definition: AMReX_Array4.H:238