Skip to content

Commit

Permalink
Implement changes by Maya Zeldin in a separate commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@RemiLehe
RemiLehe committed Sep 3, 2024
1 parent 8c4f1d4 commit 20d1fba
Show file tree
Hide file tree
Showing 7 changed files with 405 additions and 2 deletions.
146 changes: 145 additions & 1 deletion Source/Initialization/InjectorMomentum.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include <AMReX_Config.H>
#include <AMReX_Dim3.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_GpuAsyncArray.H>
#include <AMReX_REAL.H>
#include <AMReX_Parser.H>
#include <AMReX_Random.H>
Expand Down Expand Up @@ -155,6 +156,101 @@ private:
int m_flux_direction;
};

#ifdef AMREX_USE_EB
// struct whose getMomentum returns momentum for 1 particle, from random
// gaussian flux distribution in the specified direction.
// Direction is determined by stl file plane.
struct InjectorMomentumSTLGaussianFlux
{
InjectorMomentumSTLGaussianFlux (amrex::Real a_un_m, amrex::Real a_un_th,
amrex::AsyncArray<amrex::Box>& a_ba,
amrex::AsyncArray<amrex::Array4<const amrex::Real>>& a_arrays,
int a_size) noexcept
: u_m(a_un_m), u_th(a_un_th), ba(a_ba), arrays(a_arrays), size(a_size)
{
}

[[nodiscard]]
AMREX_GPU_HOST_DEVICE
amrex::XDim3
getMomentum (amrex::Real x, amrex::Real y, amrex::Real z,
amrex::RandomEngine const& engine) const noexcept
{
amrex::Real u = generateGaussianFluxDist(u_m, u_th, engine);
return calcMomentum(x, y, z, u);
}

[[nodiscard]]
AMREX_GPU_HOST_DEVICE
amrex::XDim3
getBulkMomentum (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
{
return calcMomentum(x, y, z, u_m);
}


private:

[[nodiscard]]
AMREX_GPU_HOST_DEVICE
amrex::XDim3
calcMomentum (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real momentum) const noexcept
{
#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
amrex::Box const* b_data = ba.gpuData();
amrex::Array4<const amrex::Real> const* n_data = arrays.gpuData();
return calcMomentumHelper(x, y, z, momentum, b_data, n_data);
#else
amrex::Box const* b_data = ba.data();
amrex::Array4<const amrex::Real> const* n_data = arrays.data();
return calcMomentumHelper(x, y, z, momentum, b_data, n_data);
#endif
}


[[nodiscard]]
AMREX_GPU_HOST_DEVICE
amrex::XDim3
calcMomentumHelper (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real momentum,
amrex::Box const* b_data, amrex::Array4<const amrex::Real> const* n_data) const noexcept
{

using namespace amrex::literals;

if ((b_data == nullptr) || (n_data == nullptr)) {
return amrex::XDim3{0.0_rt, 0.0_rt, 0.0_rt};
}

// how can i use x y z given that ba and array:4 have to take in integers
int int_x = amrex::Math::round(x);
int int_y = amrex::Math::round(y);
int int_z = amrex::Math::round(z);

for (int i = 0; i < size; ++i) {
if (b_data[i].contains(int_x, int_y, int_z)) {
amrex::Array4<const amrex::Real> const& eb_bnd_normal_arr = n_data[i];

amrex::Real m_i = eb_bnd_normal_arr(int_x, int_y, int_z, 1);
amrex::Real m_j = eb_bnd_normal_arr(int_x, int_y, int_z, 2);
amrex::Real m_k = eb_bnd_normal_arr(int_x, int_y, int_z, 3);

// add a variable that determines if it should go inwards or outwards (default is inwards)
amrex::Real const ux = m_i * momentum * -1;
amrex::Real const uy = m_j * momentum * -1;
amrex::Real const uz = m_k * momentum * -1;
return amrex::XDim3{ux, uy, uz};
}
}
return amrex::XDim3{0.0_rt, 0.0_rt, 0.0_rt};
}

amrex::Real u_m, u_th;
amrex::AsyncArray<amrex::Box>& ba;
amrex::AsyncArray<amrex::Array4<const amrex::Real>>& arrays;
int size;
};

#endif

// struct whose getMomentum returns momentum for 1 particle, from random
// uniform distribution u_min < u < u_max
Expand Down Expand Up @@ -540,6 +636,18 @@ struct InjectorMomentum
object(t,a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th)
{ }

#ifdef AMREX_USE_EB
// This constructor stores a InjectorMomentumSTLGaussianFlux in union object.
InjectorMomentum (InjectorMomentumSTLGaussianFlux* t,
amrex::Real a_un_m, amrex::Real a_un_th,
amrex::AsyncArray<amrex::Box>& a_ba,
amrex::AsyncArray<amrex::Array4<const amrex::Real>>& a_arrays,
int a_size)
: type(Type::stlgaussianflux),
object(t, a_un_m, a_un_th, a_ba, a_arrays, a_size)
{ }
#endif

// This constructor stores a InjectorMomentumGaussianFlux in union object.
InjectorMomentum (InjectorMomentumGaussianFlux* t,
amrex::Real a_ux_m, amrex::Real a_uy_m, amrex::Real a_uz_m,
Expand Down Expand Up @@ -615,6 +723,12 @@ struct InjectorMomentum
{
return object.gaussianflux.getMomentum(x,y,z,engine);
}
#ifdef AMREX_USE_EB
case Type::stlgaussianflux:
{
return object.stlgaussianflux.getMomentum(x,y,z,engine);
}
#endif
case Type::uniform:
{
return object.uniform.getMomentum(x,y,z,engine);
Expand Down Expand Up @@ -668,6 +782,12 @@ struct InjectorMomentum
{
return object.gaussianflux.getBulkMomentum(x,y,z);
}
#ifdef AMREX_USE_EB
case Type::stlgaussianflux:
{
return object.stlgaussianflux.getBulkMomentum(x,y,z);
}
#endif
case Type::uniform:
{
return object.uniform.getBulkMomentum(x,y,z);
Expand Down Expand Up @@ -696,7 +816,20 @@ struct InjectorMomentum
}
}

enum struct Type { constant, gaussian, gaussianflux, uniform, boltzmann, juttner, radial_expansion, parser, gaussianparser };
enum struct Type {
constant,
gaussian,
gaussianflux,
#ifdef AMREX_USE_EB
stlgaussianflux,
#endif
uniform,
boltzmann,
juttner,
radial_expansion,
parser,
gaussianparser
};
Type type;

private:
Expand All @@ -719,6 +852,14 @@ private:
amrex::Real a_uy_th, amrex::Real a_uz_th,
int a_flux_normal_axis, int a_flux_direction) noexcept
: gaussianflux(a_ux_m,a_uy_m,a_uz_m,a_ux_th,a_uy_th,a_uz_th,a_flux_normal_axis,a_flux_direction) {}
#ifdef AMREX_USE_EB
Object (InjectorMomentumSTLGaussianFlux*,
amrex::Real a_un_m, amrex::Real a_un_th,
amrex::AsyncArray<amrex::Box>& a_ba,
amrex::AsyncArray<amrex::Array4<const amrex::Real>>& a_arrays,
int a_size) noexcept
: stlgaussianflux(a_un_m, a_un_th, a_ba, a_arrays, a_size) {}
#endif
Object (InjectorMomentumUniform*,
amrex::Real a_ux_min, amrex::Real a_uy_min,
amrex::Real a_uz_min, amrex::Real a_ux_max,
Expand Down Expand Up @@ -750,6 +891,9 @@ private:
InjectorMomentumConstant constant;
InjectorMomentumGaussian gaussian;
InjectorMomentumGaussianFlux gaussianflux;
#ifdef AMREX_USE_EB
InjectorMomentumSTLGaussianFlux stlgaussianflux;
#endif
InjectorMomentumUniform uniform;
InjectorMomentumBoltzmann boltzmann;
InjectorMomentumJuttner juttner;
Expand Down
3 changes: 3 additions & 0 deletions Source/Initialization/InjectorMomentum.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,9 @@ void InjectorMomentum::clear () // NOLINT(readability-make-member-function-const
case Type::gaussian:
case Type::gaussianparser:
case Type::gaussianflux:
#ifdef AMREX_USE_EB
case Type::stlgaussianflux:
#endif
case Type::uniform:
case Type::boltzmann:
case Type::juttner:
Expand Down
106 changes: 105 additions & 1 deletion Source/Initialization/InjectorPosition.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,87 @@ private:
int dir;
};

// struct whose getPositionUnitBox returns x, y and z for a particle with
// random distribution on a plane inside a unit cell.
struct InjectorPositionRandomSTLPlane
{
InjectorPositionRandomSTLPlane (amrex::AsyncArray<amrex::Box>& a_ba,
amrex::GpuArray<amrex::Array4<const amrex::Real>, 12>& a_carray) noexcept
: ba(a_ba), carray(a_carray) {}

[[nodiscard]]
AMREX_GPU_HOST_DEVICE
amrex::XDim3
getPositionUnitBox (int /*i_part*/, amrex::IntVect const /*ref_fac*/,
amrex::RandomEngine const& engine) const noexcept
{
#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
printf("getting box in InjPosSTLPlane\n");
amrex::Box const* b_data = ba.gpuData();
printf("carray size is: %d\n", carray.size());
return getPositionHelper(b_data, engine);
#else
std::cout << "getting position in InjPosSTLPlane" << std::endl;
amrex::Box const* b_data = ba.data();
return getPositionHelper(b_data, engine);
#endif
}

private:

[[nodiscard]]
AMREX_GPU_HOST_DEVICE
amrex::XDim3
getPositionHelper (amrex::Box const* b_data,
amrex::RandomEngine const& engine) const noexcept
{

using namespace amrex::literals;

#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
printf("getting box in InjPosSTLPlane Helper\n");
#endif

// choose a random value from b_data
amrex::Box const& box = b_data[0];

#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
printf("getting eb_cent in InjPosSTLPlane Helper\n");
#endif

// get mins and maxs of the x, y, z values of the randomly selected cell given
const auto lo = amrex::lbound(box);
const auto hi = amrex::lbound(box);

int x_size = hi.x - lo.x;
int y_size = hi.y - lo.y;
int z_size = hi.z - lo.z;

#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
printf("getting random indexes in InjPosSTLPlane Helper\n");
#endif
int index_x = lo.x + amrex::Random_int(x_size, engine);
int index_y = lo.y + amrex::Random_int(y_size, engine);
int index_z = lo.z + amrex::Random_int(z_size, engine);


#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
printf("getting (x, y, z) values\n");
#endif

// getting the 0th array4 just for simplicity for now
// should just be able to get a random index
amrex::Array4<const amrex::Real> array_c = carray[0];
amrex::Real x = array_c(index_x, index_y, index_z, 1);
amrex::Real y = array_c(index_x, index_y, index_z, 2);
amrex::Real z = array_c(index_x, index_y, index_z, 3);
return amrex::XDim3{x, y, z};
}

amrex::AsyncArray<amrex::Box>& ba;
amrex::GpuArray<amrex::Array4<const amrex::Real>, 12>& carray;
};

// struct whose getPositionUnitBox returns x, y and z for a particle with
// regular distribution inside a unit cell.
struct InjectorPositionRegular
Expand Down Expand Up @@ -144,6 +225,20 @@ struct InjectorPosition
zmin(a_zmin), zmax(a_zmax)
{ }

// This constructor stores a InjectorPositionRandomSTLPlane in union object.
InjectorPosition (InjectorPositionRandomSTLPlane* t,
amrex::Real a_xmin, amrex::Real a_xmax,
amrex::Real a_ymin, amrex::Real a_ymax,
amrex::Real a_zmin, amrex::Real a_zmax,
amrex::AsyncArray<amrex::Box>& a_ba,
amrex::GpuArray<amrex::Array4<const amrex::Real>, 12>& a_carray)
: type(Type::randomstlplane),
object(t, a_ba, a_carray),
xmin(a_xmin), xmax(a_xmax),
ymin(a_ymin), ymax(a_ymax),
zmin(a_zmin), zmax(a_zmax)
{ }

// This constructor stores a InjectorPositionRegular in union object.
InjectorPosition (InjectorPositionRegular* t,
amrex::Real a_xmin, amrex::Real a_xmax,
Expand Down Expand Up @@ -180,6 +275,10 @@ struct InjectorPosition
{
return object.regular.getPositionUnitBox(i_part, ref_fac, engine);
}
case Type::randomstlplane:
{
return object.randomstlplane.getPositionUnitBox(i_part, ref_fac, engine);
}
case Type::randomplane:
{
return object.randomplane.getPositionUnitBox(i_part, ref_fac, engine);
Expand Down Expand Up @@ -233,7 +332,7 @@ struct InjectorPosition
}

private:
enum struct Type { random, randomplane, regular };
enum struct Type { random, randomplane, randomstlplane, regular };
Type type;

// An instance of union Object constructs and stores any one of
Expand All @@ -242,10 +341,15 @@ private:
Object (InjectorPositionRandom*) noexcept : random() {}
Object (InjectorPositionRandomPlane*, int const& a_dir) noexcept
: randomplane(a_dir) {}
Object (InjectorPositionRandomSTLPlane*,
amrex::AsyncArray<amrex::Box>& a_ba,
amrex::GpuArray<amrex::Array4<const amrex::Real>, 12>& a_carray) noexcept
: randomstlplane(a_ba, a_carray) {}
Object (InjectorPositionRegular*, amrex::Dim3 const& a_ppc) noexcept
: regular(a_ppc) {}
InjectorPositionRandom random;
InjectorPositionRandomPlane randomplane;
InjectorPositionRandomSTLPlane randomstlplane;
InjectorPositionRegular regular;
};
Object object;
Expand Down
1 change: 1 addition & 0 deletions Source/Initialization/PlasmaInjector.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,7 @@ protected:
void setupNuniformPerCell (amrex::ParmParse const& pp_species);
void setupExternalFile (amrex::ParmParse const& pp_species);

void setupSTLFluxInjection (amrex::ParmParse const& pp_species, const amrex::Geometry& geom);
void parseFlux (amrex::ParmParse const& pp_species);
};

Expand Down
Loading

0 comments on commit 20d1fba

Please sign in to comment.