WarpX
NodalFieldGather.H
Go to the documentation of this file.
1 /* Copyright 2019-2022 Modern Electron, Axel Huebl, Remi Lehe
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef ABLASTR_NODALFIELDGATHER_H_
8 #define ABLASTR_NODALFIELDGATHER_H_
9 
10 #include "ablastr/utils/TextMsg.H"
11 
12 #include <AMReX_Array.H>
13 #include <AMReX_Extension.H>
14 #include <AMReX_GpuQualifiers.H>
15 #include <AMReX_Math.H>
16 #include <AMReX_REAL.H>
17 
18 namespace ablastr::particles
19 {
34 void compute_weights (const amrex::ParticleReal xp,
35  const amrex::ParticleReal yp,
36  const amrex::ParticleReal zp,
39  int& i, int& j, int& k, amrex::Real W[AMREX_SPACEDIM][2], int nodal=1) noexcept
40 {
41  using namespace amrex::literals;
42 
43 #if !((nodal==0)||(nodal==1))
44  ABLASTR_ABORT_WITH_MESSAGE("Error: 'nodal' has to be equal to 0 or 1");
45 #endif
46 
47 #if (defined WARPX_DIM_3D)
48  const amrex::Real x = (xp - plo[0]) * dxi[0] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
49  const amrex::Real y = (yp - plo[1]) * dxi[1] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
50  const amrex::Real z = (zp - plo[2]) * dxi[2] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
51 
52  i = static_cast<int>(amrex::Math::floor(x));
53  j = static_cast<int>(amrex::Math::floor(y));
54  k = static_cast<int>(amrex::Math::floor(z));
55 
56  W[0][1] = x - i;
57  W[1][1] = y - j;
58  W[2][1] = z - k;
59 
60  W[0][0] = 1.0_rt - W[0][1];
61  W[1][0] = 1.0_rt - W[1][1];
62  W[2][0] = 1.0_rt - W[2][1];
63 
64 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
65 
66 # if (defined WARPX_DIM_XZ)
67  const amrex::Real x = (xp - plo[0]) * dxi[0] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
69  i = static_cast<int>(amrex::Math::floor(x));
70  W[0][1] = x - i;
71 # elif (defined WARPX_DIM_RZ)
72  const amrex::Real r = (std::sqrt(xp*xp+yp*yp) - plo[0]) * dxi[0] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
73  i = static_cast<int>(amrex::Math::floor(r));
74  W[0][1] = r - i;
75 # endif
76 
77  const amrex::Real z = (zp - plo[1]) * dxi[1] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
78  j = static_cast<int>(amrex::Math::floor(z));
79  W[1][1] = z - j;
80 
81  W[0][0] = 1.0_rt - W[0][1];
82  W[1][0] = 1.0_rt - W[1][1];
83 
84  k = 0;
85 #else
86  amrex::ignore_unused(xp, yp, zp, plo, dxi, i, j, k, W, nodal);
87  ABLASTR_ABORT_WITH_MESSAGE("Error: compute_weights not yet implemented in 1D");
88 #endif
89 }
90 
99 amrex::Real interp_field_nodal (int i, int j, int k,
100  const amrex::Real W[AMREX_SPACEDIM][2],
101  amrex::Array4<const amrex::Real> const& scalar_field) noexcept
102 {
103 #if (defined WARPX_DIM_3D)
104  amrex::Real value = 0;
105  value += scalar_field(i, j , k ) * W[0][0] * W[1][0] * W[2][0];
106  value += scalar_field(i+1, j , k ) * W[0][1] * W[1][0] * W[2][0];
107  value += scalar_field(i, j+1, k ) * W[0][0] * W[1][1] * W[2][0];
108  value += scalar_field(i+1, j+1, k ) * W[0][1] * W[1][1] * W[2][0];
109  value += scalar_field(i, j , k+1) * W[0][0] * W[1][0] * W[2][1];
110  value += scalar_field(i+1, j , k+1) * W[0][1] * W[1][0] * W[2][1];
111  value += scalar_field(i , j+1, k+1) * W[0][0] * W[1][1] * W[2][1];
112  value += scalar_field(i+1, j+1, k+1) * W[0][1] * W[1][1] * W[2][1];
113 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
114  amrex::Real value = 0;
115  value += scalar_field(i, j , k) * W[0][0] * W[1][0];
116  value += scalar_field(i+1, j , k) * W[0][1] * W[1][0];
117  value += scalar_field(i, j+1, k) * W[0][0] * W[1][1];
118  value += scalar_field(i+1, j+1, k) * W[0][1] * W[1][1];
119 #else
120  const amrex::Real value = 0;
121  amrex::ignore_unused(i, j, k, W, scalar_field);
122  ABLASTR_ABORT_WITH_MESSAGE("Error: interp_field not yet implemented in 1D");
123 #endif
124  return value;
125 }
126 
137 amrex::Real doGatherScalarFieldNodal (const amrex::ParticleReal xp,
138  const amrex::ParticleReal yp,
139  const amrex::ParticleReal zp,
140  amrex::Array4<const amrex::Real> const& scalar_field,
143 {
144  // first find the weight of surrounding nodes to use during interpolation
145  int ii, jj, kk;
146  amrex::Real W[AMREX_SPACEDIM][2];
147  compute_weights(xp, yp, zp, lo, dxi, ii, jj, kk, W);
148 
149  return interp_field_nodal(ii, jj, kk, W, scalar_field);
150 }
151 
163 doGatherVectorFieldNodal (const amrex::ParticleReal xp,
164  const amrex::ParticleReal yp,
165  const amrex::ParticleReal zp,
166  amrex::Array4<const amrex::Real> const& vector_field_x,
167  amrex::Array4<const amrex::Real> const& vector_field_y,
168  amrex::Array4<const amrex::Real> const& vector_field_z,
171 {
172  // first find the weight of surrounding nodes to use during interpolation
173  int ii, jj, kk;
174  amrex::Real W[AMREX_SPACEDIM][2];
175  compute_weights(xp, yp, zp, lo, dxi, ii, jj, kk, W);
176 
177  amrex::GpuArray<amrex::Real, 3> const field_interp = {
178  interp_field_nodal(ii, jj, kk, W, vector_field_x),
179  interp_field_nodal(ii, jj, kk, W, vector_field_y),
180  interp_field_nodal(ii, jj, kk, W, vector_field_z)
181  };
182 
183  return field_interp;
184 }
185 
186 } // namespace ablastr::particles
187 
188 #endif // ABLASTR_NODALFIELDGATHER_H_
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
#define ABLASTR_ABORT_WITH_MESSAGE(MSG)
Definition: TextMsg.H:77
Definition: DepositCharge.H:21
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::Real interp_field_nodal(int i, int j, int k, const amrex::Real W[AMREX_SPACEDIM][2], amrex::Array4< const amrex::Real > const &scalar_field) noexcept
Interpolate nodal field value based on surrounding indices and weights.
Definition: NodalFieldGather.H:99
AMREX_GPU_HOST_DEVICE AMREX_INLINE void compute_weights(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &plo, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &dxi, int &i, int &j, int &k, amrex::Real W[AMREX_SPACEDIM][2], int nodal=1) noexcept
Compute weight of each surrounding node (or cell-centered nodes) in interpolating a nodal ((or a cell...
Definition: NodalFieldGather.H:34
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::Real doGatherScalarFieldNodal(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, amrex::Array4< const amrex::Real > const &scalar_field, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &dxi, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &lo) noexcept
Scalar field gather for a single particle. The field has to be defined at the cell nodes (see https:/...
Definition: NodalFieldGather.H:137
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::GpuArray< amrex::Real, 3 > doGatherVectorFieldNodal(const amrex::ParticleReal xp, const amrex::ParticleReal yp, const amrex::ParticleReal zp, amrex::Array4< const amrex::Real > const &vector_field_x, amrex::Array4< const amrex::Real > const &vector_field_y, amrex::Array4< const amrex::Real > const &vector_field_z, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &dxi, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &lo) noexcept
Vector field gather for a single particle. The field has to be defined at the cell nodes (see https:/...
Definition: NodalFieldGather.H:163
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
i
Definition: check_interp_points_and_weights.py:174
ii
Definition: check_interp_points_and_weights.py:148
jj
Definition: check_interp_points_and_weights.py:160
value
Definition: updateAMReX.py:141