Skip to content

GK sheath BC: adding a mu dependence to parallel cut velocity #960

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Antoinehoff merged 5 commits into from
May 19, 2026
Merged

GK sheath BC: adding a mu dependence to parallel cut velocity #960

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
82efe9b
Update the kernels with `float(expand(` instead of `float(fullratsimp…
Antoinehoff May 15, 2026
6025a4e
add temp variable for vcut_fact dimensions and rename acquisition fun…
Antoinehoff May 15, 2026
2a0a198
Merge remote-tracking branch 'origin/main' into gk_sheath_vcut_mu_dep
Antoinehoff May 15, 2026
4d8fb31
The GPU implementation was done on the sheath_ai branch but not on th…
Antoinehoff May 15, 2026
b65b451
define and use vcut_fact_dim instead of cdim
Antoinehoff May 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1,556 changes: 861 additions & 695 deletions gyrokinetic/ker/bc_sheath_gyrokinetic/bc_sheath_gyrokinetic_ser_p1.c
View file
Open in desktop

Large diffs are not rendered by default.

16 changes: 8 additions & 8 deletions gyrokinetic/ker/bc_sheath_gyrokinetic/gkyl_bc_sheath_gyrokinetic_kernels.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,14 @@ static inline double invL(double x) {

EXTERN_C_BEG

GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_lower_1x1v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_upper_1x1v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_lower_1x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_upper_1x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_lower_2x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_upper_2x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_lower_3x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_upper_3x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_lower_1x1v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_upper_1x1v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_lower_1x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_upper_1x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_lower_2x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_upper_2x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_lower_3x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);
GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_upper_3x2v_ser_p1(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);

EXTERN_C_END

305 changes: 219 additions & 86 deletions gyrokinetic/unit/ctest_bc_sheath_gyrokinetic.c
View file
Open in desktop

Large diffs are not rendered by default.

74 changes: 70 additions & 4 deletions gyrokinetic/zero/bc_sheath_gyrokinetic.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,17 @@
#include <gkyl_alloc.h>
#include <assert.h>

static struct gkyl_array*
mkarr(bool on_gpu, long nc, long size)
{
struct gkyl_array* a;
if (on_gpu)
a = gkyl_array_cu_dev_new(GKYL_DOUBLE, nc, size);
else
a = gkyl_array_new(GKYL_DOUBLE, nc, size);
return a;
}

struct gkyl_bc_sheath_gyrokinetic*
gkyl_bc_sheath_gyrokinetic_new(int dir, enum gkyl_edge_loc edge, const struct gkyl_basis *basis,
const struct gkyl_range *skin_r, const struct gkyl_range *ghost_r, const struct gkyl_velocity_map *vel_map,
Expand All @@ -21,9 +32,36 @@ gkyl_bc_sheath_gyrokinetic_new(int dir, enum gkyl_edge_loc edge, const struct gk
up->skin_r = skin_r;
up->ghost_r = ghost_r;
up->vel_map = gkyl_velocity_map_acquire(vel_map);
int vdim = skin_r->ndim - cdim;

// Choose the kernel that does the reflection/no reflection/partial
// reflection.
up->vcut_fact_dim = cdim-1 + vdim-1;

// Get polynomial order
int poly_order;
if (use_gpu) {
struct gkyl_basis *basis_ho;
basis_ho = gkyl_malloc(sizeof(struct gkyl_basis));
gkyl_cu_memcpy(basis_ho, basis, sizeof(struct gkyl_basis), GKYL_CU_MEMCPY_D2H);
poly_order = basis_ho->poly_order;
gkyl_free(basis_ho);
} else {
poly_order = basis->poly_order;
}
gkyl_cart_modal_serendip(&up->vcut_fact_basis, up->vcut_fact_dim, poly_order);

int lower[up->vcut_fact_dim], upper[up->vcut_fact_dim];
for (int d=0; d < cdim-1; d++) {
lower[d] = skin_r->lower[d];
upper[d] = skin_r->upper[d];
}
lower[up->vcut_fact_dim-1] = skin_r->lower[skin_r->ndim-1];
upper[up->vcut_fact_dim-1] = skin_r->upper[skin_r->ndim-1];
gkyl_range_init(&up->vcut_fact_local, up->vcut_fact_basis.ndim, lower, upper);
up->vcut_fact = mkarr(use_gpu, up->vcut_fact_basis.num_basis, up->vcut_fact_local.volume);
double dg_norm = pow(sqrt(2), up->vcut_fact_basis.ndim);
gkyl_array_shiftc_range(up->vcut_fact, dg_norm, 0, &up->vcut_fact_local);

// Choose the kernels that does the reflection/no reflection/partial reflection,
up->kernels = gkyl_malloc(sizeof(struct gkyl_bc_sheath_gyrokinetic_kernels));
#ifdef GKYL_HAVE_CUDA
if (use_gpu) {
Expand Down Expand Up @@ -57,6 +95,7 @@ gkyl_bc_sheath_gyrokinetic_advance(const struct gkyl_bc_sheath_gyrokinetic *up,

int fidx[GKYL_MAX_DIM]; // Flipped index.
int vidx[2];
int vcut_fact_idx[up->vcut_fact_local.ndim];

int pdim = up->skin_r->ndim;
int vpar_dir = up->cdim;
Expand All @@ -80,24 +119,50 @@ gkyl_bc_sheath_gyrokinetic_advance(const struct gkyl_bc_sheath_gyrokinetic *up,
for (int d=up->cdim; d<pdim; d++) vidx[d-up->cdim] = iter.idx[d];
long conf_loc = gkyl_range_idx(conf_r, iter.idx);
long vel_loc = gkyl_range_idx(&up->vel_map->local_vel, vidx);

const double *phi_p = (const double*) gkyl_array_cfetch(phi, conf_loc);
const double *phi_wall_p = (const double*) gkyl_array_cfetch(phi_wall, conf_loc);
const double *vmap_p = (const double*) gkyl_array_cfetch(up->vel_map->vmap, vel_loc);

int vcut_fact_dim = up->vcut_fact_local.ndim;
for (int d=0; d<vcut_fact_dim-1; d++) vcut_fact_idx[d] = iter.idx[d]; // config space perp directions.
vcut_fact_idx[vcut_fact_dim-1] = iter.idx[pdim-1]; // mu direction.
long vcut_fact_loc = gkyl_range_idx(&up->vcut_fact_local, vcut_fact_idx);
const double *vcut_fact_p = (const double*) gkyl_array_cfetch(up->vcut_fact, vcut_fact_loc);

// Calculate reflected distribution function fhat.
// note: reflected distribution can be
// 1) fhat=0 (no reflection, i.e. absorb),
// 2) fhat=f (full reflection)
// 3) fhat=c*f (partial reflection)
double fhat[up->basis->num_basis];
up->kernels->reflectedf(vmap_p, up->q2Dm, phi_p, phi_wall_p, inp, fhat);
up->kernels->reflectedf(vmap_p, up->q2Dm, phi_p, phi_wall_p, vcut_fact_p, inp, fhat);

// Reflect fhat into skin cells.
bc_gksheath_reflect(up->dir, up->basis, up->cdim, out, fhat);
}
}

void gkyl_bc_sheath_gyrokinetic_set_vcut_fact(const struct gkyl_bc_sheath_gyrokinetic *up, const struct gkyl_array *vcut_fact)
{
gkyl_array_copy_range(up->vcut_fact, vcut_fact, &up->vcut_fact_local);
}

struct gkyl_array* gkyl_bc_sheath_gyrokinetic_acquire_vcut_fact(struct gkyl_bc_sheath_gyrokinetic *up)
{
return gkyl_array_acquire(up->vcut_fact);
}

struct gkyl_basis gkyl_bc_sheath_gyrokinetic_get_vcut_fact_basis(struct gkyl_bc_sheath_gyrokinetic *up)
{
return up->vcut_fact_basis;
}

struct gkyl_range* gkyl_bc_sheath_gyrokinetic_get_vcut_fact_range(struct gkyl_bc_sheath_gyrokinetic *up)
{
return &up->vcut_fact_local;
}

void gkyl_bc_sheath_gyrokinetic_release(struct gkyl_bc_sheath_gyrokinetic *up)
{
// Release memory associated with this updater.
Expand All @@ -107,6 +172,7 @@ void gkyl_bc_sheath_gyrokinetic_release(struct gkyl_bc_sheath_gyrokinetic *up)
}
#endif
gkyl_velocity_map_release(up->vel_map);
gkyl_array_release(up->vcut_fact);
gkyl_free(up->kernels);
gkyl_free(up);
}
23 changes: 17 additions & 6 deletions gyrokinetic/zero/bc_sheath_gyrokinetic_cu.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,14 +32,17 @@ gkyl_bc_gksheath_choose_reflectedf_kernel_cu(const struct gkyl_basis *basis,
}

__global__ static void
gkyl_bc_sheath_gyrokinetic_advance_cu_ker(int cdim, int dir, const struct gkyl_range skin_r, const struct gkyl_range ghost_r,
const struct gkyl_range conf_r, const struct gkyl_range vel_r, const struct gkyl_basis *basis,
const struct gkyl_array *vmap, double q2Dm, const struct gkyl_array *phi,
const struct gkyl_array *phi_wall, struct gkyl_bc_sheath_gyrokinetic_kernels *kers, struct gkyl_array *distf)
gkyl_bc_sheath_gyrokinetic_advance_cu_ker(int cdim, int dir, const struct gkyl_range skin_r,
const struct gkyl_range ghost_r, const struct gkyl_range conf_r, const struct gkyl_range vcut_r,
const struct gkyl_range vel_r, const struct gkyl_basis *basis, const struct gkyl_array *vmap,
double q2Dm, const struct gkyl_array *phi, const struct gkyl_array *phi_wall,
const struct gkyl_array *vcut_fact, struct gkyl_bc_sheath_gyrokinetic_kernels *kers,
struct gkyl_array *distf)
{
int fidx[GKYL_MAX_DIM]; // Flipped index.
int pidx[GKYL_MAX_DIM];
int vidx[2];
int vcut_fact_idx[GKYL_MAX_CDIM];

int pdim = skin_r.ndim;
int vpar_dir = cdim;
Expand All @@ -53,6 +56,7 @@ gkyl_bc_sheath_gyrokinetic_advance_cu_ker(int cdim, int dir, const struct gkyl_r
// since update_range is a subrange

gkyl_sub_range_inv_idx(&skin_r, linc, pidx);
// gkyl_sub_range_inv_idx(&vcut_r, linc, vcut_fact_idx);

gkyl_copy_int_arr(pdim, pidx, fidx);
fidx[vpar_dir] = uplo - pidx[vpar_dir];
Expand All @@ -73,13 +77,20 @@ gkyl_bc_sheath_gyrokinetic_advance_cu_ker(int cdim, int dir, const struct gkyl_r
const double *phi_wall_p = (const double*) gkyl_array_cfetch(phi_wall, conf_loc);
const double *vmap_p = (const double*) gkyl_array_cfetch(vmap, vel_loc);

// Get vcut factor.
int vcut_fact_dim = vcut_r.ndim;
for (int d=0; d<vcut_fact_dim-1; d++) vcut_fact_idx[d] = pidx[d]; // config space perp directions.
vcut_fact_idx[vcut_fact_dim-1] = pidx[pdim-1]; // mu direction.
long vcut_fact_loc = gkyl_range_idx(&vcut_r, vcut_fact_idx);
const double *vcut_fact_p = (const double*) gkyl_array_cfetch(vcut_fact, vcut_fact_loc);

// Calculate reflected distribution function fhat.
// note: reflected distribution can be
// 1) fhat=0 (no reflection, i.e. absorb),
// 2) fhat=f (full reflection)
// 3) fhat=c*f (partial reflection)
double fhat[112]; // MF 2022/08/24: hardcoded to number of DG coeffs in 3x2v p2 for now.
kers->reflectedf(vmap_p, q2Dm, phi_p, phi_wall_p, inp, fhat);
kers->reflectedf(vmap_p, q2Dm, phi_p, phi_wall_p, vcut_fact_p, inp, fhat);

// Reflect fhat into skin cells.
bc_gksheath_reflect(dir, basis, cdim, out, fhat);
Expand All @@ -94,7 +105,7 @@ gkyl_bc_sheath_gyrokinetic_advance_cu(const struct gkyl_bc_sheath_gyrokinetic *u
int nblocks = up->skin_r->nblocks, nthreads = up->skin_r->nthreads;

gkyl_bc_sheath_gyrokinetic_advance_cu_ker<<<nblocks, nthreads>>>(up->cdim, up->dir, *up->skin_r, *up->ghost_r,
*conf_r, up->vel_map->local_vel, up->basis, up->vel_map->vmap->on_dev, up->q2Dm, phi->on_dev, phi_wall->on_dev,
*conf_r, up->vcut_fact_local, up->vel_map->local_vel, up->basis, up->vel_map->vmap->on_dev, up->q2Dm, phi->on_dev, phi_wall->on_dev, up->vcut_fact->on_dev,
up->kernels_cu, distf->on_dev);
}
}
35 changes: 35 additions & 0 deletions gyrokinetic/zero/gkyl_bc_sheath_gyrokinetic.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,12 +33,47 @@ struct gkyl_bc_sheath_gyrokinetic* gkyl_bc_sheath_gyrokinetic_new(int dir, enum
* @param up BC updater.
* @param phi Electrostatic potential.
* @param phi_wall Wall potential.
* @param vcut_fact Alpha parameter for sheath BCs.
* @param distf Distribution function array to apply BC to.
* @param conf_r Configuration space range (to index phi).
*/
void gkyl_bc_sheath_gyrokinetic_advance(const struct gkyl_bc_sheath_gyrokinetic *up, const struct gkyl_array *phi,
const struct gkyl_array *phi_wall, struct gkyl_array *distf, const struct gkyl_range *conf_r);

/**
* Set the vcut_fact array used in the sheath BC.
* This is used to implement a mu-dependent vcut in the sheath BC,
* where vcut_fact is the mu dependent multiplying factor in the expression for vcut.
*
* @param up BC updater.
* @param vcut_fact The vcut_fact array to use in the sheath BC.
*/
void gkyl_bc_sheath_gyrokinetic_set_vcut_fact(const struct gkyl_bc_sheath_gyrokinetic *up, const struct gkyl_array *vcut_fact);

/**
* Acquire the vcut_fact array used in the sheath BC. The pointer needs to be released by the caller using gkyl_array_release.
*
* @param up BC updater.
* @return The vcut_fact array used in the sheath BC.
*/
struct gkyl_array* gkyl_bc_sheath_gyrokinetic_acquire_vcut_fact(struct gkyl_bc_sheath_gyrokinetic *up);

/**
* Get the basis of the vcut_fact array used in the sheath BC.
*
* @param up BC updater.
* @return The basis of the vcut_fact array used in the sheath BC.
*/
struct gkyl_basis gkyl_bc_sheath_gyrokinetic_get_vcut_fact_basis(struct gkyl_bc_sheath_gyrokinetic *up);

/**
* Get the range of the vcut_fact array used in the sheath BC.
*
* @param up BC updater.
* @return The range of the vcut_fact array used in the sheath BC.
*/
struct gkyl_range* gkyl_bc_sheath_gyrokinetic_get_vcut_fact_range(struct gkyl_bc_sheath_gyrokinetic *up);

/**
* Free memory associated with bc_sheath_gyrokinetic updater.
*
Expand Down
6 changes: 5 additions & 1 deletion gyrokinetic/zero/gkyl_bc_sheath_gyrokinetic_priv.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@

// Function pointer type for sheath reflection kernels.
typedef void (*sheath_reflectedf_t)(const double *vmap, const double q2Dm,
const double *phi, const double *phiWall, const double *f, double *fRefl);
const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl);

typedef struct { sheath_reflectedf_t kernels[3]; } sheath_reflectedf_kern_list; // For use in kernel tables.
typedef struct { sheath_reflectedf_kern_list list[4]; } edged_sheath_reflectedf_kern_list;
Expand Down Expand Up @@ -48,6 +48,10 @@ struct gkyl_bc_sheath_gyrokinetic {
struct gkyl_bc_sheath_gyrokinetic_kernels *kernels_cu; // device copy.
const struct gkyl_range *skin_r, *ghost_r; // Skin and ghost ranges.
const struct gkyl_velocity_map *vel_map; // Velocity space mapping.
int vcut_fact_dim; // Dimensionality of vcut_fact array.
struct gkyl_array *vcut_fact; // factor for mu dependent vcut in sheath BC.
struct gkyl_basis vcut_fact_basis; // Basis for vcut_fact array (expansion in perpendicular config space and mu).
struct gkyl_range vcut_fact_local; // Range for vcut_fact array.
};

void
Expand Down
Loading