From 281e77d458af2bc7270a959b9d7772cdc791a91c Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Thu, 2 Jul 2026 11:49:20 -0400 Subject: [PATCH 01/12] Review and change a bit how the FLR effects are handled. A factor 1/2 was missing in the species FLR operators. The gk_field_flr is also changed and now does nothing as it should invert the flr on phi. There is now a placeholder where we will call the application for the operator 1 -rho^2 laplacian_perp which is currently not possible. --- core/zero/gkyl_eqn_type.h | 3 +- gyrokinetic/apps/gk_field_1x.c | 3 + gyrokinetic/apps/gk_field_flr.c | 84 ++++++++++-------------- gyrokinetic/apps/gk_species.c | 6 +- gyrokinetic/apps/gkyl_gyrokinetic.h | 4 ++ gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 6 +- 6 files changed, 50 insertions(+), 56 deletions(-) diff --git a/core/zero/gkyl_eqn_type.h b/core/zero/gkyl_eqn_type.h index b3388eab11..cd345c3366 100644 --- a/core/zero/gkyl_eqn_type.h +++ b/core/zero/gkyl_eqn_type.h @@ -168,7 +168,8 @@ enum gkyl_vel_edge { // Identifiers for FLR models (in gyrokinetics). enum gkyl_gk_flr_type { GKYL_GK_FLR_NONE = 0, // No FLR effects. - GKYL_GK_FLR_PADE_CONST, // Pade-based approx. w/ const. rho_ts=sqrt(Tperp_s/m_s) + GKYL_GK_FLR_PADE_CONST, // Pade-based approx. w/ const. rho_ts; at the field level, uses a single reference gyroradius in the operator retrieving phi. + GKYL_GK_FLR_PADE_CONST_SUM, // Pade-based approx. w/ const. rho_ts; at the field level, uses the polarization-weighted average of the species gyroradii. }; // Gyrokinetic anomaous diffusion models. diff --git a/gyrokinetic/apps/gk_field_1x.c b/gyrokinetic/apps/gk_field_1x.c index 2d91e55a4e..81b85b0160 100644 --- a/gyrokinetic/apps/gk_field_1x.c +++ b/gyrokinetic/apps/gk_field_1x.c @@ -17,6 +17,9 @@ gk_field_rhs_phi_1x(struct gkyl_gyrokinetic_app *app, struct gk_field *field) { // Solve the Poisson equation in 1x with the parallel FEM projection. gk_field_fem_projection_par(app, field, field->rho_c, field->phi_smooth); + + // Finish the Poisson solve with FLR effects. + field->invert_flr(app, field, field->phi_smooth); } static void diff --git a/gyrokinetic/apps/gk_field_flr.c b/gyrokinetic/apps/gk_field_flr.c index 10c068fb05..2cbc164db5 100644 --- a/gyrokinetic/apps/gk_field_flr.c +++ b/gyrokinetic/apps/gk_field_flr.c @@ -12,59 +12,49 @@ gk_field_flr_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f) assert(app->cdim > 1); f->invert_flr = gk_field_invert_flr; - double flr_weight = 0.0; - for (int i = 0; i < app->num_species; ++i){ - struct gk_species *s = &app->species[i]; - double gyroradius_bmag = s->info.flr.bmag ? s->info.flr.bmag : app->bmag_ref; - flr_weight += s->info.flr.Tperp * s->info.mass / (pow(s->info.charge * gyroradius_bmag, 2.0)); + // Reference (squared) gyroradius in the operator A = 1 - rho^2*nabla_perp^2 + // used to retrieve phi from the modified potential Phi_0 (step 4 of the + // algorithm in flr_effects.tex). + double rhoSq_ref = 0.0; + if (f->info.flr.type == GKYL_GK_FLR_PADE_CONST) { + // Single reference species (e.g. the main ion). + assert(f->info.flr.gyroradius > 0.0); + rhoSq_ref = pow(f->info.flr.gyroradius, 2.0); } - // Initialize the weight in the Laplacian operator. - f->flr_rhoSq_sum = mkarr(app->use_gpu, (2 * (app->cdim - 1) - 1) * app->basis.num_basis, app->local_ext.volume); - gkyl_array_set_offset(f->flr_rhoSq_sum, flr_weight, app->gk_geom->geo_int.gxxj, 0 * app->basis.num_basis); - if (app->cdim > 2) { - gkyl_array_set_offset(f->flr_rhoSq_sum, flr_weight, app->gk_geom->geo_int.gxyj, 1 * app->basis.num_basis); - gkyl_array_set_offset(f->flr_rhoSq_sum, flr_weight, app->gk_geom->geo_int.gyyj, 2 * app->basis.num_basis); - } - // Initialize the factor multiplying the field in the FLR operator. - f->flr_kSq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); - gkyl_array_shiftc(f->flr_kSq, -pow(sqrt(2.0), app->cdim), 0); // Sets kSq=-1. - - // If domain is not periodic use Dirichlet BCs. - struct gkyl_poisson_bc flr_bc = {}; - - bool bc_is_np[GKYL_MAX_CDIM]; // Is the BC in this direction non-periodic? - for (int d = 0; d < app->cdim; ++d) { - bc_is_np[d] = true; - } - for (int d = 0; d < app->num_periodic_dir; ++d) { - bc_is_np[app->periodic_dirs[d]] = false; - } - - for (int d = 0; d < app->cdim - 1; d++) { - if (bc_is_np[d]) { - struct gkyl_gyrokinetic_bc *bc_lo = gk_fetch_bc_with_dir_edge(f->info.poisson_bcs, 2 * app->cdim, d, GKYL_LOWER_EDGE); - if (bc_lo != 0) { - flr_bc.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_DIRICHLET_VARYING); - } - - struct gkyl_gyrokinetic_bc *bc_up = gk_fetch_bc_with_dir_edge(f->info.poisson_bcs, 2 * app->cdim, d, GKYL_UPPER_EDGE); - if (bc_up != 0) { - flr_bc.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_DIRICHLET_VARYING); - } - } else { - flr_bc.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); - flr_bc.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); + else { + // GKYL_GK_FLR_PADE_CONST_SUM, or field-level FLR options not set (default): + // polarization-weighted average of the species gyroradii, + // rho^2 = sum_s eps_s0*rho_s0^2 / sum_s eps_s0, eps_s0 = n_s0*m_s/B^2, + // which suppresses the electron contribution by m_e/m_i and reduces to + // rho_i0^2 for a single ion species. Species without FLR enabled + // contribute rho_s0^2 = 0. + double polarization_bmag = f->info.polarization_bmag ? f->info.polarization_bmag : app->bmag_ref; + double eps_sum = 0.0; + for (int i = 0; i < app->num_species; ++i) { + struct gk_species *s = &app->species[i]; + double gyroradius_bmag = s->info.flr.bmag ? s->info.flr.bmag : app->bmag_ref; + double rhoSq_s = s->info.flr.Tperp * s->info.mass / pow(s->info.charge * gyroradius_bmag, 2.0); + double eps_s0 = s->info.polarization_density * s->info.mass / pow(polarization_bmag, 2.0); + rhoSq_ref += eps_s0 * rhoSq_s; + eps_sum += eps_s0; } + rhoSq_ref /= eps_sum; + } + // Weight in the perpendicular Laplacian of A, staged for the forthcoming + // Laplacian-apply infrastructure. + f->flr_rhoSq = mkarr(app->use_gpu, (2 * (app->cdim - 1) - 1) * app->basis.num_basis, app->local_ext.volume); + gkyl_array_set_offset(f->flr_rhoSq, rhoSq_ref, app->gk_geom->geo_int.gxxj, 0 * app->basis.num_basis); + if (app->cdim > 2) { + gkyl_array_set_offset(f->flr_rhoSq, rhoSq_ref, app->gk_geom->geo_int.gxyj, 1 * app->basis.num_basis); + gkyl_array_set_offset(f->flr_rhoSq, rhoSq_ref, app->gk_geom->geo_int.gyyj, 2 * app->basis.num_basis); } - // Deflated Poisson solve is performed on range assuming decomposition is *only* in z. - f->flr_op = gkyl_deflated_fem_poisson_new(app->grid, app->basis_on_dev, app->basis, - app->local, app->local, f->flr_rhoSq_sum, f->flr_kSq, flr_bc, NULL, app->use_gpu); } void gk_field_invert_flr(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array *phi) { - gkyl_deflated_fem_poisson_advance(field->flr_op, phi, phi, phi); + // Retrieve phi from the modified potential Phi_0 by applying + // A = 1 - rho_i0^2*nabla_perp^2 } void @@ -75,7 +65,5 @@ gk_field_invert_flr_none(gkyl_gyrokinetic_app *app, struct gk_field *field, stru void gk_field_flr_release(const struct gkyl_gyrokinetic_app *app, struct gk_field *f) { - gkyl_array_release(f->flr_rhoSq_sum); - gkyl_array_release(f->flr_kSq); - gkyl_deflated_fem_poisson_release(f->flr_op); + gkyl_array_release(f->flr_rhoSq); } diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index 2b93be63cf..4e54695719 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/gyrokinetic/apps/gk_species.c @@ -1564,10 +1564,10 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st double flr_weight = gks->info.flr.Tperp*gks->info.mass/(pow(gks->info.charge*gyroradius_bmag,2.0)); // Initialize the weight in the Laplacian operator. gks->flr_rhoSqD2 = mkarr(app->use_gpu, (2*(app->cdim-1)-1)*app->basis.num_basis, app->local_ext.volume); - gkyl_array_set_offset(gks->flr_rhoSqD2, flr_weight, app->gk_geom->geo_int.gxxj, 0*app->basis.num_basis); + gkyl_array_set_offset(gks->flr_rhoSqD2, 0.5*flr_weight, app->gk_geom->geo_int.gxxj, 0*app->basis.num_basis); if (app->cdim > 2) { - gkyl_array_set_offset(gks->flr_rhoSqD2, flr_weight, app->gk_geom->geo_int.gxyj, 1*app->basis.num_basis); - gkyl_array_set_offset(gks->flr_rhoSqD2, flr_weight, app->gk_geom->geo_int.gyyj, 2*app->basis.num_basis); + gkyl_array_set_offset(gks->flr_rhoSqD2, 0.5*flr_weight, app->gk_geom->geo_int.gxyj, 1*app->basis.num_basis); + gkyl_array_set_offset(gks->flr_rhoSqD2, 0.5*flr_weight, app->gk_geom->geo_int.gyyj, 2*app->basis.num_basis); } // Initialize the factor multiplying the field in the FLR operator. gks->flr_kSq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h index e04e1510ba..0b07cf3dc8 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic.h @@ -335,6 +335,8 @@ struct gkyl_gyrokinetic_flr { double Tperp; // Perp temperature used to evaluate gyroradius. double bmag; // Magnetic field used to evaluate gyroradius. If not provided // it'll use B in the center of the domain. + double gyroradius; // Reference (e.g. main ion) gyroradius used in the field-level + // operator retrieving phi from the modified potential. }; struct gkyl_gyrokinetic_correct_inp { @@ -583,6 +585,8 @@ struct gkyl_gyrokinetic_field { bool phi_wall_up_evolve; // set to true if biased wall potential on upper wall function is time dependent struct gkyl_poisson_bias_line_list *bias_line_list; // Biased lines constraining the solution. + + struct gkyl_gyrokinetic_flr flr; // Options for FLR effects. }; struct gkyl_gyrokinetic_eirene { diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index 5615160973..6fb8bf7a48 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -1354,10 +1354,8 @@ struct gk_field { // Objects needed for FLR effects. bool use_flr; // Whether to apply FLR effects. - void (*invert_flr)(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array *phi); // Function inverting FLR operator. - struct gkyl_array *flr_rhoSq_sum; // Laplacian weight in FLR operator. - struct gkyl_array *flr_kSq; // Field multiplying phi in FLR operator. - struct gkyl_deflated_fem_poisson *flr_op; // Helmholtz solver to invert FLR operator. + void (*invert_flr)(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array *phi); // Function retrieving phi from the modified potential Phi_0 by inverting the FLR operator, i.e. applying A = 1 - rho^2*nabla_perp^2. + struct gkyl_array *flr_rhoSq; // rho^2 weight (times J*g^ij) in the perpendicular Laplacian of A. struct gkyl_array_integrate *calc_em_energy; // Operator computing EM energy. double *em_energy_red, *em_energy_red_global; // memory for use in GPU reduction of EM energy From 0677f1b94c06a2bc0d6fe57558840b560353295b Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Thu, 2 Jul 2026 18:32:44 -0400 Subject: [PATCH 02/12] Implementation of a way to apply the LHS of a screened poisson problem to a DG array. Basically, while the current poisson fem solver solves (M+K) phi = M rho for phi, the new routine gkyl_fem_poisson_lhs_apply computes rho = M^-1 (M+K) phi which is the DG representation of the linear operator applied to phi. This is done using superlu sp_dgemv routine on CPU. There is no equivalent in cuDSS so we use cusparse cusparse SpMV routine on GPU. Cusparse provides a routine that allows to points on an existing matrix, cusparseCreateCsr, which prevents double allocation/computation of the matrix in the poisson updater. The unit test associated to this new feature tests that we can revert a poisson solve by applying the operator back. I need to test the GPU code on perlmutter. --- core/zero/cudss_ops.cu | 34 ++++ core/zero/cusolver_ops.cu | 34 ++++ core/zero/gkyl_cudss_ops.h | 10 ++ core/zero/gkyl_cusolver_ops.h | 17 ++ core/zero/gkyl_superlu_ops.h | 10 ++ core/zero/superlu_ops.c | 7 + moments/unit/ctest_fem_poisson_lhs_apply.c | 200 +++++++++++++++++++++ moments/zero/fem_poisson.c | 72 ++++++++ moments/zero/fem_poisson_cu.cu | 16 ++ moments/zero/gkyl_fem_poisson.h | 13 ++ moments/zero/gkyl_fem_poisson_priv.h | 12 +- 11 files changed, 424 insertions(+), 1 deletion(-) create mode 100644 moments/unit/ctest_fem_poisson_lhs_apply.c diff --git a/core/zero/cudss_ops.cu b/core/zero/cudss_ops.cu index 47284f3711..40382e7069 100644 --- a/core/zero/cudss_ops.cu +++ b/core/zero/cudss_ops.cu @@ -1,6 +1,7 @@ #ifdef GKYL_HAVE_CUDSS #include +#include extern "C" { #include @@ -263,6 +264,39 @@ gkyl_culinsolver_sync(struct gkyl_culinsolver_prob *prob) cudaStreamSynchronize(prob->stream); } +void +gkyl_culinsolver_mat_vec(struct gkyl_culinsolver_prob *prob, const double *x, double *y) +{ + cusparseHandle_t handle; + cusparseCreate(&handle); + cusparseSetStream(handle, prob->stream); + + cusparseSpMatDescr_t matA; + cusparseCreateCsr(&matA, prob->mrow, prob->ncol, prob->nnz, + prob->csr_rowptr_cu, prob->csr_colind_cu, prob->csr_val_cu, + CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F); + + cusparseDnVecDescr_t vecX, vecY; + cusparseCreateDnVec(&vecX, prob->ncol, (void*) x, CUDA_R_64F); + cusparseCreateDnVec(&vecY, prob->mrow, y, CUDA_R_64F); + + double alpha = 1.0, beta = 0.0; + size_t buf_sz = 0; + cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, + CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &buf_sz); + void *buf = buf_sz > 0? gkyl_cu_malloc(buf_sz) : NULL; + + cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, + CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, buf); + cudaStreamSynchronize(prob->stream); + + if (buf) gkyl_cu_free(buf); + cusparseDestroyDnVec(vecX); + cusparseDestroyDnVec(vecY); + cusparseDestroySpMat(matA); + cusparseDestroy(handle); +} + void gkyl_culinsolver_finish_host(struct gkyl_culinsolver_prob *prob) { diff --git a/core/zero/cusolver_ops.cu b/core/zero/cusolver_ops.cu index 8634a30ed8..d9449d2705 100644 --- a/core/zero/cusolver_ops.cu +++ b/core/zero/cusolver_ops.cu @@ -425,6 +425,40 @@ gkyl_culinsolver_solve(struct gkyl_culinsolver_prob *prob) cusolverRfBatchSolve(prob->cusolverRfH, prob->d_P, prob->d_Q, 1, prob->d_T, prob->mrow, prob->rhspointers_cu, prob->mrow); } +void +gkyl_culinsolver_sync(struct gkyl_culinsolver_prob *prob) +{ + cudaStreamSynchronize(prob->stream); +} + +void +gkyl_culinsolver_mat_vec(struct gkyl_culinsolver_prob *prob, const double *x, double *y) +{ + cusparseSpMatDescr_t matA; + cusparseCreateCsr(&matA, prob->mrow, prob->ncol, prob->nnz, + prob->csrrowptrA_cu, prob->csrcolindA_cu, prob->csrvalA_cu, + CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F); + + cusparseDnVecDescr_t vecX, vecY; + cusparseCreateDnVec(&vecX, prob->ncol, (void*) x, CUDA_R_64F); + cusparseCreateDnVec(&vecY, prob->mrow, y, CUDA_R_64F); + + double alpha = 1.0, beta = 0.0; + size_t buf_sz = 0; + cusparseSpMV_bufferSize(prob->cusparseH, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, + CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &buf_sz); + void *buf = buf_sz > 0? gkyl_cu_malloc(buf_sz) : NULL; + + cusparseSpMV(prob->cusparseH, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, + CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, buf); + cudaStreamSynchronize(prob->stream); + + if (buf) gkyl_cu_free(buf); + cusparseDestroyDnVec(vecX); + cusparseDestroyDnVec(vecY); + cusparseDestroySpMat(matA); +} + void gkyl_culinsolver_finish_host(struct gkyl_culinsolver_prob *prob) { diff --git a/core/zero/gkyl_cudss_ops.h b/core/zero/gkyl_cudss_ops.h index 239457be02..90ad0dc9ee 100644 --- a/core/zero/gkyl_cudss_ops.h +++ b/core/zero/gkyl_cudss_ops.h @@ -66,6 +66,16 @@ void gkyl_culinsolver_brhs_from_triples(struct gkyl_culinsolver_prob *prob, gkyl */ void gkyl_culinsolver_solve(struct gkyl_culinsolver_prob *prob); +/** + * Compute the matrix-vector product y = A*x (on the device) for the (single) + * problem matrix A. Assumes nprob=1. Uses cuSPARSE on the CSR arrays. + * + * @param prob cuDSS struct holding the assembled A matrix. + * @param x Input vector on the device (length ncol). + * @param y Output vector on the device (length mrow). + */ +void gkyl_culinsolver_mat_vec(struct gkyl_culinsolver_prob *prob, const double *x, double *y); + /** * Copy solution back to host * diff --git a/core/zero/gkyl_cusolver_ops.h b/core/zero/gkyl_cusolver_ops.h index ade556c3e9..86a1fa103c 100644 --- a/core/zero/gkyl_cusolver_ops.h +++ b/core/zero/gkyl_cusolver_ops.h @@ -51,6 +51,23 @@ void gkyl_culinsolver_brhs_from_triples(struct gkyl_culinsolver_prob *prob, gkyl */ void gkyl_culinsolver_solve(struct gkyl_culinsolver_prob *prob); +/** + * Synchronize the stream the solver ran on. + * + * @param prob cuSolver struct holding arrays used in problem. + */ +void gkyl_culinsolver_sync(struct gkyl_culinsolver_prob *prob); + +/** + * Compute the matrix-vector product y = A*x (on the device) for the (single) + * problem matrix A. Assumes nprob=1. Uses cuSPARSE on the CSR arrays. + * + * @param prob cuSolver struct holding the assembled A matrix. + * @param x Input vector on the device (length ncol). + * @param y Output vector on the device (length mrow). + */ +void gkyl_culinsolver_mat_vec(struct gkyl_culinsolver_prob *prob, const double *x, double *y); + /** * Copy solution back to host * diff --git a/core/zero/gkyl_superlu_ops.h b/core/zero/gkyl_superlu_ops.h index c69cbb2f97..5faa10fdd6 100644 --- a/core/zero/gkyl_superlu_ops.h +++ b/core/zero/gkyl_superlu_ops.h @@ -65,6 +65,16 @@ void gkyl_superlu_brhs_from_array(struct gkyl_superlu_prob *prob, const double * */ void gkyl_superlu_solve(struct gkyl_superlu_prob *prob); +/** + * Compute the matrix-vector product y = A*x for the (single) problem matrix A. + * Assumes nprob=1. Does not require A to be factorized. + * + * @param prob SuperLu struct holding the assembled A matrix. + * @param x Input vector (length ncol). + * @param y Output vector (length mrow). + */ +void gkyl_superlu_mat_vec(struct gkyl_superlu_prob *prob, const double *x, double *y); + /** * Update the values of the left-hand-side SuperLU matrix A in Ax=B problem from a list of * triples. Assumes the sparsity pattern didn't change. diff --git a/core/zero/superlu_ops.c b/core/zero/superlu_ops.c index ffbfe2094d..612b7169e6 100644 --- a/core/zero/superlu_ops.c +++ b/core/zero/superlu_ops.c @@ -319,6 +319,13 @@ gkyl_superlu_solve(struct gkyl_superlu_prob *prob) } } +void +gkyl_superlu_mat_vec(struct gkyl_superlu_prob *prob, const double *x, double *y) +{ + char trans[2] = "N"; + sp_dgemv(trans, 1.0, prob->A[0], (double *) x, 1, 0.0, y, 1); +} + void gkyl_superlu_amat_update_from_triples(struct gkyl_superlu_prob *prob, struct gkyl_mat_triples **tri) { diff --git a/moments/unit/ctest_fem_poisson_lhs_apply.c b/moments/unit/ctest_fem_poisson_lhs_apply.c new file mode 100644 index 0000000000..bacae2c84a --- /dev/null +++ b/moments/unit/ctest_fem_poisson_lhs_apply.c @@ -0,0 +1,200 @@ +// Test gkyl_fem_poisson_lhs_apply, which applies the LHS operator +// M^{-1}*(M+K) = weak (1 - rho^2*Lap) using the assembled FEM matrices. +// +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +static struct gkyl_array* +mkarr(long nc, long size) +{ + return gkyl_array_new(GKYL_DOUBLE, nc, size); +} + +static struct gkyl_array* +mkarr_dev(bool on_gpu, long nc, long size) +{ +#ifdef GKYL_HAVE_CUDA + if (on_gpu) return gkyl_array_cu_dev_new(GKYL_DOUBLE, nc, size); +#endif + return gkyl_array_new(GKYL_DOUBLE, nc, size); +} + +static double +calc_l2(struct gkyl_rect_grid grid, struct gkyl_range range, struct gkyl_range range_ext, + struct gkyl_basis basis, struct gkyl_array *f1, struct gkyl_array *f2) +{ + struct gkyl_array *diff = mkarr(basis.num_basis, range_ext.volume); + struct gkyl_range_iter iter; + gkyl_range_iter_init(&iter, &range); + while (gkyl_range_iter_next(&iter)) { + long lidx = gkyl_range_idx(&range, iter.idx); + const double *f1p = gkyl_array_cfetch(f1, lidx); + const double *f2p = f2? gkyl_array_cfetch(f2, lidx) : NULL; + double *diffp = gkyl_array_fetch(diff, lidx); + for (int i=0; i 2) { + gkyl_array_shiftc(epsilon_ho, rho*rho*dg0norm, 2*basis.num_basis); + } + struct gkyl_array *kSq_ho = mkarr(basis.num_basis, local_ext.volume); + gkyl_array_shiftc(kSq_ho, -dg0norm, 0); + + struct gkyl_array *src_ho = mkarr(basis.num_basis, local_ext.volume); + gkyl_proj_on_basis *proj = gkyl_proj_on_basis_new(&grid, &basis, poly_order+1, 1, + dim==1? evalFunc_1x : evalFunc_2x, NULL); + gkyl_proj_on_basis_advance(proj, 0.0, &local, src_ho); + gkyl_proj_on_basis_release(proj); + + // Working arrays (on device if use_gpu), populated from the host copies. + struct gkyl_array *epsilon = mkarr_dev(use_gpu, neps*basis.num_basis, local_ext.volume); + struct gkyl_array *kSq = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + struct gkyl_array *src = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + struct gkyl_array *w = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + struct gkyl_array *g = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + struct gkyl_array *w2 = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + gkyl_array_copy(epsilon, epsilon_ho); + gkyl_array_copy(kSq, kSq_ho); + gkyl_array_copy(src, src_ho); + + gkyl_fem_poisson *poisson = gkyl_fem_poisson_new(&local, &grid, basis, &bcs, NULL, + epsilon, kSq, false, use_gpu); + + // Continuous field w = A^{-1}*src. + gkyl_fem_poisson_set_rhs(poisson, src, NULL); + gkyl_fem_poisson_solve(poisson, w); + + // Apply the LHS operator, g = (1 - rho^2*Lap) w. + gkyl_fem_poisson_lhs_apply(poisson, w, g); + + // Exact identity: solving with the same operator inverts the apply. + gkyl_fem_poisson_set_rhs(poisson, g, NULL); + gkyl_fem_poisson_solve(poisson, w2); + + // Copy results to host for the L2 checks. + struct gkyl_array *w_ho = mkarr(basis.num_basis, local_ext.volume); + struct gkyl_array *g_ho = mkarr(basis.num_basis, local_ext.volume); + struct gkyl_array *w2_ho = mkarr(basis.num_basis, local_ext.volume); + gkyl_array_copy(w_ho, w); + gkyl_array_copy(g_ho, g); + gkyl_array_copy(w2_ho, w2); + + double err = calc_l2(grid, local, local_ext, basis, w2_ho, w_ho); + double norm_w = calc_l2(grid, local, local_ext, basis, w_ho, NULL); + TEST_CHECK(err < 1.0e-10*norm_w); + TEST_MSG("round-trip L2 error = %g (|w| = %g)", err, norm_w); + + // Applying (1 - rho^2*Lap) amplifies the mode, so |g| > |w|. + double norm_g = calc_l2(grid, local, local_ext, basis, g_ho, NULL); + TEST_CHECK(norm_g > norm_w); + TEST_MSG("|g| = %g, |w| = %g", norm_g, norm_w); + + gkyl_array_release(epsilon_ho); + gkyl_array_release(kSq_ho); + gkyl_array_release(src_ho); + gkyl_array_release(epsilon); + gkyl_array_release(kSq); + gkyl_array_release(src); + gkyl_array_release(w); + gkyl_array_release(g); + gkyl_array_release(w2); + gkyl_array_release(w_ho); + gkyl_array_release(g_ho); + gkyl_array_release(w2_ho); + gkyl_fem_poisson_release(poisson); +} + +void test_1x_p1(void) { + int cells[] = {16}; + test_lhs_apply(1, cells, 0.3, GKYL_POISSON_DIRICHLET, false); + test_lhs_apply(1, cells, 0.3, GKYL_POISSON_PERIODIC, false); +} + +void test_2x_p1(void) { int cells[] = {24, 16}; + test_lhs_apply(2, cells, 0.3, GKYL_POISSON_DIRICHLET, false); + test_lhs_apply(2, cells, 0.3, GKYL_POISSON_PERIODIC, false); + } + +#ifdef GKYL_HAVE_CUDA + +void gpu_test_1x_p1(void) { + int cells[] = {16}; + test_lhs_apply(1, cells, 0.3, GKYL_POISSON_DIRICHLET, true); + test_lhs_apply(1, cells, 0.3, GKYL_POISSON_PERIODIC, true); +} + +void gpu_test_2x_p1(void) { + int cells[] = {24, 16}; + test_lhs_apply(2, cells, 0.3, GKYL_POISSON_DIRICHLET, true); + test_lhs_apply(2, cells, 0.3, GKYL_POISSON_PERIODIC, true); +} + +#endif + +TEST_LIST = { + { "test_1x_p1", test_1x_p1 }, + { "test_2x_p1", test_2x_p1 }, +#ifdef GKYL_HAVE_CUDA + { "gpu_test_1x_p1", gpu_test_1x_p1 }, + { "gpu_test_2x_p1", gpu_test_2x_p1 }, +#endif + { NULL, NULL }, +}; diff --git a/moments/zero/fem_poisson.c b/moments/zero/fem_poisson.c index 43682a56c3..01941a72c6 100644 --- a/moments/zero/fem_poisson.c +++ b/moments/zero/fem_poisson.c @@ -1,6 +1,7 @@ #include #include #include +#include static void fem_poisson_bias_src_disabled(gkyl_fem_poisson* up, struct gkyl_array *rhsin) @@ -61,6 +62,8 @@ gkyl_fem_poisson_new(const struct gkyl_range *solve_range, const struct gkyl_rec up->kernels_cu = up->kernels; #endif + up->has_lhs_apply = false; + up->bcs = *bcs; up->solve_range = solve_range; up->ndim = grid->ndim; up->grid = *grid; @@ -403,8 +406,77 @@ gkyl_fem_poisson_solve(gkyl_fem_poisson* up, struct gkyl_array *phiout) { } +static void +fem_poisson_lhs_apply_init(gkyl_fem_poisson *up) +{ + // Build a mass-matrix (M) solver: same operator with epsilon=0 and a constant + // kSq=-1 (so the LHS stencil reduces to the mass matrix). + struct gkyl_array *eps_zero, *kSq_mass; +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) { + eps_zero = gkyl_array_cu_dev_new(GKYL_DOUBLE, up->epsilon->ncomp, up->epsilon->size); + kSq_mass = gkyl_array_cu_dev_new(GKYL_DOUBLE, up->num_basis, up->epsilon->size); + } else { + eps_zero = gkyl_array_new(GKYL_DOUBLE, up->epsilon->ncomp, up->epsilon->size); + kSq_mass = gkyl_array_new(GKYL_DOUBLE, up->num_basis, up->epsilon->size); + } +#else + eps_zero = gkyl_array_new(GKYL_DOUBLE, up->epsilon->ncomp, up->epsilon->size); + kSq_mass = gkyl_array_new(GKYL_DOUBLE, up->num_basis, up->epsilon->size); +#endif + gkyl_array_clear(eps_zero, 0.0); + gkyl_array_clear(kSq_mass, 0.0); + gkyl_array_shiftc(kSq_mass, -pow(sqrt(2.0), up->ndim), 0); + + up->mass = gkyl_fem_poisson_new(up->solve_range, &up->grid, up->basis, &up->bcs, NULL, + eps_zero, kSq_mass, true, up->use_gpu); + + up->lhs_dual = gkyl_malloc(sizeof(double[up->numnodes_global])); +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) + up->lhs_dual_cu = gkyl_cu_malloc(sizeof(double[up->numnodes_global])); +#endif + + gkyl_array_release(eps_zero); + gkyl_array_release(kSq_mass); + up->has_lhs_apply = true; +} + +void +gkyl_fem_poisson_lhs_apply(gkyl_fem_poisson* up, struct gkyl_array *xin, struct gkyl_array *xout) +{ + assert(up->ishelmholtz); + if (!up->has_lhs_apply) fem_poisson_lhs_apply_init(up); + +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) { + gkyl_fem_poisson_lhs_apply_cu(up, xin, xout); + return; + } +#endif + + // Recover the global nodal vector x_nodal = M^{-1}*(M_src*xin). + gkyl_fem_poisson_set_rhs(up->mass, xin, NULL); + gkyl_superlu_solve(up->mass->prob); + double *x_nodal = gkyl_superlu_get_rhs_ptr(up->mass->prob, 0); + + // Dual d = (M+K)*x_nodal, then solve M*z = d (-> modal xout). + gkyl_superlu_mat_vec(up->prob, x_nodal, up->lhs_dual); + gkyl_superlu_brhs_from_array(up->mass->prob, up->lhs_dual); + gkyl_fem_poisson_solve(up->mass, xout); +} + void gkyl_fem_poisson_release(gkyl_fem_poisson *up) { + if (up->has_lhs_apply) { + gkyl_fem_poisson_release(up->mass); + gkyl_free(up->lhs_dual); +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) + gkyl_cu_free(up->lhs_dual_cu); +#endif + } + if (up->isdomperiodic) { gkyl_array_release(up->rhs_cellavg); gkyl_free(up->rhs_avg); diff --git a/moments/zero/fem_poisson_cu.cu b/moments/zero/fem_poisson_cu.cu index 41f69050c7..5c893740ac 100644 --- a/moments/zero/fem_poisson_cu.cu +++ b/moments/zero/fem_poisson_cu.cu @@ -319,3 +319,19 @@ gkyl_fem_poisson_solve_cu(gkyl_fem_poisson *up, struct gkyl_array *phiout) x_cu, *up->solve_range, up->kernels_cu); } +void +gkyl_fem_poisson_lhs_apply_cu(gkyl_fem_poisson *up, struct gkyl_array *xin, struct gkyl_array *xout) +{ + // Recover the global nodal vector x_nodal = M^{-1}*(M_src*xin). + gkyl_fem_poisson_set_rhs(up->mass, xin, NULL); + gkyl_culinsolver_solve(up->mass->prob_cu); + gkyl_culinsolver_sync(up->mass->prob_cu); + double *x_nodal = gkyl_culinsolver_get_sol_ptr(up->mass->prob_cu, 0); + + // Dual d = (M+K)*x_nodal, then solve M*z = d (-> modal xout). + gkyl_culinsolver_mat_vec(up->prob_cu, x_nodal, up->lhs_dual_cu); + double *rhs_cu = gkyl_culinsolver_get_rhs_ptr(up->mass->prob_cu, 0); + gkyl_cu_memcpy(rhs_cu, up->lhs_dual_cu, sizeof(double)*up->numnodes_global, GKYL_CU_MEMCPY_D2D); + gkyl_fem_poisson_solve(up->mass, xout); +} + diff --git a/moments/zero/gkyl_fem_poisson.h b/moments/zero/gkyl_fem_poisson.h index 99c4847f60..96766abc03 100644 --- a/moments/zero/gkyl_fem_poisson.h +++ b/moments/zero/gkyl_fem_poisson.h @@ -58,6 +58,19 @@ void gkyl_fem_poisson_set_rhs(gkyl_fem_poisson* up, struct gkyl_array *rhsin, co */ void gkyl_fem_poisson_solve(gkyl_fem_poisson* up, struct gkyl_array *phiout); +/** + * Apply the LHS operator, i.e. compute xout = M^{-1}*(M+K)*xin where M is the + * mass matrix and M+K is the (Helmholtz) LHS matrix. With epsilon=rho^2 and + * kSq=-1 this returns the weak image of (1 - rho^2*Lap)*xin. Assumes xin is + * continuous (in the FEM space) and requires a Helmholtz solver (kSq!=NULL). + * Runs on the CPU (SuperLU) or GPU (cuDSS/cuSPARSE). + * + * @param up FEM poisson updater to run. + * @param xin Continuous DG field to apply the operator to. + * @param xout DG field holding the result. + */ +void gkyl_fem_poisson_lhs_apply(gkyl_fem_poisson* up, struct gkyl_array *xin, struct gkyl_array *xout); + /** * Delete updater. * diff --git a/moments/zero/gkyl_fem_poisson_priv.h b/moments/zero/gkyl_fem_poisson_priv.h index cfa1704492..176aca24f5 100644 --- a/moments/zero/gkyl_fem_poisson_priv.h +++ b/moments/zero/gkyl_fem_poisson_priv.h @@ -1959,13 +1959,23 @@ struct gkyl_fem_poisson { int num_bias_plane; // Number of biased planes. struct gkyl_poisson_bias_plane *bias_planes; // Biased planes. - bias_src_func_t bias_plane_src; // Function to enforce biasing in RHS source. + bias_src_func_t bias_plane_src; // Function to enforce biasing in RHS source. + + // Objects to apply the LHS operator (M+K), used to apply 1-rho^2*Lap. + bool has_lhs_apply; // Whether the LHS-apply objects are allocated. + struct gkyl_poisson_bc bcs; // Saved BCs (used to build the mass solver). + struct gkyl_fem_poisson *mass; // Mass-matrix (M) solver for the LHS apply. + double *lhs_dual; // Host global dual (M+K)*x vector. + double *lhs_dual_cu; // Device global dual (M+K)*x vector. }; void fem_poisson_choose_kernels_cu(const struct gkyl_basis* basis, const struct gkyl_poisson_bc* bcs, bool isvareps, const bool *isdirperiodic, struct gkyl_fem_poisson_kernels *kers); +void +gkyl_fem_poisson_lhs_apply_cu(gkyl_fem_poisson *up, struct gkyl_array *xin, struct gkyl_array *xout); + static long gkyl_fem_poisson_global_num_nodes(const int dim, const int poly_order, const int basis_type, const int *num_cells, bool *isdirperiodic) From f26260885751095c6ce81ec8d37d7d205bab01e2 Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Mon, 6 Jul 2026 09:57:31 -0400 Subject: [PATCH 03/12] Implement the apply LHS to the poisson perp operator too. The former commit was using the global fem_poisson operator which is not coherent with how we solve the field equation in GK. A new apply mat operator is created to apply matrices with the gk_fem_poisson_perp now. An equivalent unit test is provided. We can see if we want to remove the other apply LHS later. We also move poisson_bcs inside the field structure so that we can reuse it for FLR effects. It also simplifies gk_species.c as it does not have to compute it again, the only caveat is that field must be initialized before species but I think other parts of the code assume that. --- core/zero/cudss_ops.cu | 41 ++-- core/zero/cusolver_ops.cu | 41 ++-- core/zero/gkyl_cudss_ops.h | 10 +- core/zero/gkyl_cusolver_ops.h | 10 +- core/zero/gkyl_superlu_ops.h | 10 +- core/zero/superlu_ops.c | 3 +- gyrokinetic/apps/gk_field_2x3x.c | 19 +- gyrokinetic/apps/gk_field_flr.c | 31 ++- gyrokinetic/apps/gk_species.c | 33 +-- gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 5 +- .../unit/ctest_fem_poisson_perp_lhs_apply.c | 204 ++++++++++++++++++ gyrokinetic/zero/fem_poisson_perp.c | 84 ++++++++ gyrokinetic/zero/gkyl_fem_poisson_perp.h | 13 ++ gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h | 7 + 14 files changed, 417 insertions(+), 94 deletions(-) create mode 100644 gyrokinetic/unit/ctest_fem_poisson_perp_lhs_apply.c diff --git a/core/zero/cudss_ops.cu b/core/zero/cudss_ops.cu index 40382e7069..1c8b0f84d6 100644 --- a/core/zero/cudss_ops.cu +++ b/core/zero/cudss_ops.cu @@ -271,29 +271,36 @@ gkyl_culinsolver_mat_vec(struct gkyl_culinsolver_prob *prob, const double *x, do cusparseCreate(&handle); cusparseSetStream(handle, prob->stream); - cusparseSpMatDescr_t matA; - cusparseCreateCsr(&matA, prob->mrow, prob->ncol, prob->nnz, - prob->csr_rowptr_cu, prob->csr_colind_cu, prob->csr_val_cu, - CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F); - - cusparseDnVecDescr_t vecX, vecY; - cusparseCreateDnVec(&vecX, prob->ncol, (void*) x, CUDA_R_64F); - cusparseCreateDnVec(&vecY, prob->mrow, y, CUDA_R_64F); - double alpha = 1.0, beta = 0.0; + void *buf = NULL; size_t buf_sz = 0; - cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, - CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &buf_sz); - void *buf = buf_sz > 0? gkyl_cu_malloc(buf_sz) : NULL; + for (int k=0; knprob; k++) { + cusparseSpMatDescr_t matA; + cusparseCreateCsr(&matA, prob->mrow, prob->ncol, prob->nnz, + prob->csr_rowptr_cu, prob->csr_colind_cu, prob->csr_val_cu+k*prob->nnz, + CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F); + + cusparseDnVecDescr_t vecX, vecY; + cusparseCreateDnVec(&vecX, prob->ncol, (void*) &x[k*prob->ncol], CUDA_R_64F); + cusparseCreateDnVec(&vecY, prob->mrow, &y[k*prob->mrow], CUDA_R_64F); + + if (k == 0) { + // All problems have the same sparsity pattern, so one buffer suffices. + cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, + CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &buf_sz); + buf = buf_sz > 0? gkyl_cu_malloc(buf_sz) : NULL; + } - cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, - CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, buf); + cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, + CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, buf); + + cusparseDestroyDnVec(vecX); + cusparseDestroyDnVec(vecY); + cusparseDestroySpMat(matA); + } cudaStreamSynchronize(prob->stream); if (buf) gkyl_cu_free(buf); - cusparseDestroyDnVec(vecX); - cusparseDestroyDnVec(vecY); - cusparseDestroySpMat(matA); cusparseDestroy(handle); } diff --git a/core/zero/cusolver_ops.cu b/core/zero/cusolver_ops.cu index d9449d2705..bc431e7451 100644 --- a/core/zero/cusolver_ops.cu +++ b/core/zero/cusolver_ops.cu @@ -434,29 +434,36 @@ gkyl_culinsolver_sync(struct gkyl_culinsolver_prob *prob) void gkyl_culinsolver_mat_vec(struct gkyl_culinsolver_prob *prob, const double *x, double *y) { - cusparseSpMatDescr_t matA; - cusparseCreateCsr(&matA, prob->mrow, prob->ncol, prob->nnz, - prob->csrrowptrA_cu, prob->csrcolindA_cu, prob->csrvalA_cu, - CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F); - - cusparseDnVecDescr_t vecX, vecY; - cusparseCreateDnVec(&vecX, prob->ncol, (void*) x, CUDA_R_64F); - cusparseCreateDnVec(&vecY, prob->mrow, y, CUDA_R_64F); - double alpha = 1.0, beta = 0.0; + void *buf = NULL; size_t buf_sz = 0; - cusparseSpMV_bufferSize(prob->cusparseH, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, - CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &buf_sz); - void *buf = buf_sz > 0? gkyl_cu_malloc(buf_sz) : NULL; + for (int k=0; knprob; k++) { + cusparseSpMatDescr_t matA; + cusparseCreateCsr(&matA, prob->mrow, prob->ncol, prob->nnz, + prob->csrrowptrA_cu, prob->csrcolindA_cu, prob->csrvalA_cu+k*prob->nnz, + CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F); + + cusparseDnVecDescr_t vecX, vecY; + cusparseCreateDnVec(&vecX, prob->ncol, (void*) &x[k*prob->ncol], CUDA_R_64F); + cusparseCreateDnVec(&vecY, prob->mrow, &y[k*prob->mrow], CUDA_R_64F); + + if (k == 0) { + // All problems have the same sparsity pattern, so one buffer suffices. + cusparseSpMV_bufferSize(prob->cusparseH, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, + CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &buf_sz); + buf = buf_sz > 0? gkyl_cu_malloc(buf_sz) : NULL; + } - cusparseSpMV(prob->cusparseH, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, - CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, buf); + cusparseSpMV(prob->cusparseH, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, + CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, buf); + + cusparseDestroyDnVec(vecX); + cusparseDestroyDnVec(vecY); + cusparseDestroySpMat(matA); + } cudaStreamSynchronize(prob->stream); if (buf) gkyl_cu_free(buf); - cusparseDestroyDnVec(vecX); - cusparseDestroyDnVec(vecY); - cusparseDestroySpMat(matA); } void diff --git a/core/zero/gkyl_cudss_ops.h b/core/zero/gkyl_cudss_ops.h index 90ad0dc9ee..cf2b760d5a 100644 --- a/core/zero/gkyl_cudss_ops.h +++ b/core/zero/gkyl_cudss_ops.h @@ -67,12 +67,12 @@ void gkyl_culinsolver_brhs_from_triples(struct gkyl_culinsolver_prob *prob, gkyl void gkyl_culinsolver_solve(struct gkyl_culinsolver_prob *prob); /** - * Compute the matrix-vector product y = A*x (on the device) for the (single) - * problem matrix A. Assumes nprob=1. Uses cuSPARSE on the CSR arrays. + * Compute the matrix-vector products y_i = A_i*x_i (on the device) for each + * of the nprob problem matrices A_i. Uses cuSPARSE on the CSR arrays. * - * @param prob cuDSS struct holding the assembled A matrix. - * @param x Input vector on the device (length ncol). - * @param y Output vector on the device (length mrow). + * @param prob cuDSS struct holding the assembled A matrices. + * @param x Input vector on the device (length nprob*ncol). + * @param y Output vector on the device (length nprob*mrow). */ void gkyl_culinsolver_mat_vec(struct gkyl_culinsolver_prob *prob, const double *x, double *y); diff --git a/core/zero/gkyl_cusolver_ops.h b/core/zero/gkyl_cusolver_ops.h index 86a1fa103c..158bf62a3f 100644 --- a/core/zero/gkyl_cusolver_ops.h +++ b/core/zero/gkyl_cusolver_ops.h @@ -59,12 +59,12 @@ void gkyl_culinsolver_solve(struct gkyl_culinsolver_prob *prob); void gkyl_culinsolver_sync(struct gkyl_culinsolver_prob *prob); /** - * Compute the matrix-vector product y = A*x (on the device) for the (single) - * problem matrix A. Assumes nprob=1. Uses cuSPARSE on the CSR arrays. + * Compute the matrix-vector products y_i = A_i*x_i (on the device) for each + * of the nprob problem matrices A_i. Uses cuSPARSE on the CSR arrays. * - * @param prob cuSolver struct holding the assembled A matrix. - * @param x Input vector on the device (length ncol). - * @param y Output vector on the device (length mrow). + * @param prob cuSolver struct holding the assembled A matrices. + * @param x Input vector on the device (length nprob*ncol). + * @param y Output vector on the device (length nprob*mrow). */ void gkyl_culinsolver_mat_vec(struct gkyl_culinsolver_prob *prob, const double *x, double *y); diff --git a/core/zero/gkyl_superlu_ops.h b/core/zero/gkyl_superlu_ops.h index 5faa10fdd6..5b07570d5a 100644 --- a/core/zero/gkyl_superlu_ops.h +++ b/core/zero/gkyl_superlu_ops.h @@ -66,12 +66,12 @@ void gkyl_superlu_brhs_from_array(struct gkyl_superlu_prob *prob, const double * void gkyl_superlu_solve(struct gkyl_superlu_prob *prob); /** - * Compute the matrix-vector product y = A*x for the (single) problem matrix A. - * Assumes nprob=1. Does not require A to be factorized. + * Compute the matrix-vector products y_i = A_i*x_i for each of the nprob + * problem matrices A_i. Does not require the A_i to be factorized. * - * @param prob SuperLu struct holding the assembled A matrix. - * @param x Input vector (length ncol). - * @param y Output vector (length mrow). + * @param prob SuperLu struct holding the assembled A matrices. + * @param x Input vector (length nprob*ncol). + * @param y Output vector (length nprob*mrow). */ void gkyl_superlu_mat_vec(struct gkyl_superlu_prob *prob, const double *x, double *y); diff --git a/core/zero/superlu_ops.c b/core/zero/superlu_ops.c index 612b7169e6..e264e0d8c5 100644 --- a/core/zero/superlu_ops.c +++ b/core/zero/superlu_ops.c @@ -323,7 +323,8 @@ void gkyl_superlu_mat_vec(struct gkyl_superlu_prob *prob, const double *x, double *y) { char trans[2] = "N"; - sp_dgemv(trans, 1.0, prob->A[0], (double *) x, 1, 0.0, y, 1); + for (size_t k=0; knprob; k++) + sp_dgemv(trans, 1.0, prob->A[k], (double *) &x[k*prob->ncol], 1, 0.0, &y[k*prob->mrow], 1); } void diff --git a/gyrokinetic/apps/gk_field_2x3x.c b/gyrokinetic/apps/gk_field_2x3x.c index f64b5c7f15..f100c32b0e 100644 --- a/gyrokinetic/apps/gk_field_2x3x.c +++ b/gyrokinetic/apps/gk_field_2x3x.c @@ -532,33 +532,32 @@ gk_field_fem_new_2x3x(struct gkyl_gyrokinetic_app *app, struct gk_field *f) } // Translate input file BCs into Poisson BCs. - struct gkyl_poisson_bc poisson_bcs = { }; for (int d=0; dcdim-1; d++) { if (bc_is_np[d]) { struct gkyl_gyrokinetic_bc *bc_lo = gk_fetch_bc_with_dir_edge(f->info.poisson_bcs, 2*app->cdim, d, GKYL_LOWER_EDGE); if (bc_lo != 0) { - poisson_bcs.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(bc_lo->type); + f->poisson_bcs.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(bc_lo->type); for (int i=0; i<3; i++) { - poisson_bcs.lo_value[d].v[i] = bc_lo->value[i]; + f->poisson_bcs.lo_value[d].v[i] = bc_lo->value[i]; } } struct gkyl_gyrokinetic_bc *bc_up = gk_fetch_bc_with_dir_edge(f->info.poisson_bcs, 2*app->cdim, d, GKYL_UPPER_EDGE); if (bc_up != 0) { - poisson_bcs.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(bc_up->type); + f->poisson_bcs.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(bc_up->type); for (int i=0; i<3; i++) { - poisson_bcs.up_value[d].v[i] = bc_up->value[i]; + f->poisson_bcs.up_value[d].v[i] = bc_up->value[i]; } } } else { - poisson_bcs.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); - poisson_bcs.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); + f->poisson_bcs.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); + f->poisson_bcs.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); } } // Initialize the Poisson solver. f->fem_poisson_perp = gkyl_fem_poisson_perp_new(&app->local, &app->grid, app->basis, - &poisson_bcs, f->info.bias_line_list, f->epsilon, NULL, app->use_gpu); + &f->poisson_bcs, f->info.bias_line_list, f->epsilon, NULL, app->use_gpu); f->phi_bc = 0; f->is_dirichletvar = false; @@ -642,11 +641,11 @@ gk_field_fem_new_2x3x(struct gkyl_gyrokinetic_app *app, struct gk_field *f) if (app->gk_geom->has_LCFS) { // Updaters to enforce twist-and-shift and sheath BCs. - gk_field_2x3x_add_IWL_updaters(app, f, &poisson_bcs); + gk_field_2x3x_add_IWL_updaters(app, f, &f->poisson_bcs); } else if (f->bc_par_phi == GKYL_BC_GK_FIELD_TWISTSHIFT) { // Updaters to enforce twist-and-shift BCs. - gk_field_2x3x_add_TS_updaters(app, f, &poisson_bcs); + gk_field_2x3x_add_TS_updaters(app, f, &f->poisson_bcs); } // Set the pointer to the function that computes phi. diff --git a/gyrokinetic/apps/gk_field_flr.c b/gyrokinetic/apps/gk_field_flr.c index 2cbc164db5..09f24a6fd0 100644 --- a/gyrokinetic/apps/gk_field_flr.c +++ b/gyrokinetic/apps/gk_field_flr.c @@ -40,21 +40,32 @@ gk_field_flr_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f) } rhoSq_ref /= eps_sum; } - // Weight in the perpendicular Laplacian of A, staged for the forthcoming - // Laplacian-apply infrastructure. - f->flr_rhoSq = mkarr(app->use_gpu, (2 * (app->cdim - 1) - 1) * app->basis.num_basis, app->local_ext.volume); - gkyl_array_set_offset(f->flr_rhoSq, rhoSq_ref, app->gk_geom->geo_int.gxxj, 0 * app->basis.num_basis); - if (app->cdim > 2) { - gkyl_array_set_offset(f->flr_rhoSq, rhoSq_ref, app->gk_geom->geo_int.gxyj, 1 * app->basis.num_basis); - gkyl_array_set_offset(f->flr_rhoSq, rhoSq_ref, app->gk_geom->geo_int.gyyj, 2 * app->basis.num_basis); + // The apply passes boundary values through Dirichlet rows, so spatially + // varying Dirichlet BCs are not supported with FLR effects. + for (int d = 0; d < app->cdim - 1; d++) { + assert(f->poisson_bcs.lo_type[d] != GKYL_POISSON_DIRICHLET_VARYING); + assert(f->poisson_bcs.up_type[d] != GKYL_POISSON_DIRICHLET_VARYING); } + + // Weight in the perpendicular Laplacian of A = 1 - rho^2*nabla_perp^2. + f->flr_rhoSq = mkarr(app->use_gpu, (2*(app->cdim/3)+1)*app->basis.num_basis, app->local_ext.volume); + struct gkyl_array *Jgij[3] = {app->gk_geom->geo_int.gxxj, app->gk_geom->geo_int.gxyj, app->gk_geom->geo_int.gyyj}; + for (int i=0; icdim-2/app->cdim; i++) { + gkyl_array_set_offset(f->flr_rhoSq, rhoSq_ref, Jgij[i], i*app->basis.num_basis); + } + f->flr_kSq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); + gkyl_array_shiftc(f->flr_kSq, -pow(sqrt(2.0),app->cdim), 0); // Sets kSq=-1. + + f->flr_op = gkyl_fem_poisson_perp_new(&app->local, &app->grid, app->basis, &f->poisson_bcs, f->info.bias_line_list, f->flr_rhoSq, f->flr_kSq, app->use_gpu); + } void gk_field_invert_flr(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array *phi) { - // Retrieve phi from the modified potential Phi_0 by applying - // A = 1 - rho_i0^2*nabla_perp^2 + // Retrieve phi from the potential obtained from the FLR charge. + // phi = (1 - rho^2*nabla_perp^2) Phi_0 + gkyl_fem_poisson_perp_lhs_apply(field->flr_op, phi, phi); } void @@ -66,4 +77,6 @@ void gk_field_flr_release(const struct gkyl_gyrokinetic_app *app, struct gk_field *f) { gkyl_array_release(f->flr_rhoSq); + gkyl_array_release(f->flr_kSq); + gkyl_fem_poisson_perp_release(f->flr_op); } diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index 4e54695719..fe3acd9905 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/gyrokinetic/apps/gk_species.c @@ -29,9 +29,12 @@ void gk_species_gyroaverage_enabled(gkyl_gyrokinetic_app *app, struct gk_species *species, struct gkyl_array *field_in, struct gkyl_array *field_gyroavg) { + // Solve the gyroaveraging problem: (1 - rho_s^2 nabla_perp^2) phi_g = phi_in to + // get the gyroaveraged potential used in the Hamiltonian. struct timespec wst = gkyl_wall_clock(); // Gyroaverage input field. - gkyl_deflated_fem_poisson_advance(species->flr_op, field_in, field_in, field_gyroavg); + gkyl_fem_poisson_perp_set_rhs(species->flr_op, field_in); + gkyl_fem_poisson_perp_solve(species->flr_op, field_gyroavg); app->stat.species_gyroavg_tm += gkyl_time_diff_now_sec(wst); } @@ -1573,28 +1576,10 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st gks->flr_kSq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); gkyl_array_shiftc(gks->flr_kSq, -pow(sqrt(2.0),app->cdim), 0); // Sets kSq=-1. - // If domain is not periodic use Dirichlet BCs. - struct gkyl_poisson_bc flr_bc; - for (int d=0; dcdim-1; d++) { - struct gkyl_gyrokinetic_bc *bc_lo = gk_fetch_bc_with_dir_edge(app->field->info.poisson_bcs, 2*app->cdim, d, GKYL_LOWER_EDGE); - if (bc_lo != 0) { - if (bc_lo->type == GKYL_BC_GK_FIELD_PERIODIC) - flr_bc.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); - else - flr_bc.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_DIRICHLET_VARYING); - } - - struct gkyl_gyrokinetic_bc *bc_up = gk_fetch_bc_with_dir_edge(app->field->info.poisson_bcs, 2*app->cdim, d, GKYL_UPPER_EDGE); - if (bc_up != 0) { - if (bc_up->type == GKYL_BC_GK_FIELD_PERIODIC) - flr_bc.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); - else - flr_bc.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_DIRICHLET_VARYING); - } - } - // Deflated Poisson solve is performed on range assuming decomposition is *only* in z. - gks->flr_op = gkyl_deflated_fem_poisson_new(app->grid, app->basis_on_dev, app->basis, - app->local, app->local, gks->flr_rhoSqD2, gks->flr_kSq, flr_bc, NULL, app->use_gpu); + // Use the same solver and BCs as the field solve so the gyroaverage and + // the FLR inversion of phi use the same discretization. + gks->flr_op = gkyl_fem_poisson_perp_new(&app->local, &app->grid, app->basis, + &app->field->poisson_bcs, 0, gks->flr_rhoSqD2, gks->flr_kSq, app->use_gpu); } else { gks->gyroaverage = gk_species_gyroaverage_disabled; @@ -2060,7 +2045,7 @@ gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *gks if (gks->info.flr.type) { gkyl_array_release(gks->flr_rhoSqD2); gkyl_array_release(gks->flr_kSq); - gkyl_deflated_fem_poisson_release(gks->flr_op); + gkyl_fem_poisson_perp_release(gks->flr_op); } gks->release_func(app, gks); diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index 6fb8bf7a48..a300d749a4 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -1148,7 +1148,7 @@ struct gk_species { struct gkyl_array *m0_gyroavg; // Gyroaveraged particle density. struct gkyl_array *flr_rhoSqD2; // Laplacian weight in FLR operator. struct gkyl_array *flr_kSq; // Field multiplying phi in FLR operator. - struct gkyl_deflated_fem_poisson *flr_op; // Helmholtz solver to invert FLR operator. + struct gkyl_fem_poisson_perp *flr_op; // Screened Poisson solver used to gyroaverage fields. // Pointer to function that performs the gyroaverage. void (*gyroaverage)(gkyl_gyrokinetic_app *app, struct gk_species *species, struct gkyl_array *field_in, struct gkyl_array *field_gyroavg); @@ -1351,11 +1351,14 @@ struct gk_field { struct gkyl_fem_parproj *fem_parproj_rho_core, *fem_parproj_phi_core; // FEM projection in the core. struct gkyl_fem_poisson_perp *fem_poisson_perp; // Solves - nabla . (epsilon * nabla phi) - kSq * phi = rho. + struct gkyl_poisson_bc poisson_bcs; // Boundary conditions for Poisson solver. // Objects needed for FLR effects. bool use_flr; // Whether to apply FLR effects. void (*invert_flr)(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array *phi); // Function retrieving phi from the modified potential Phi_0 by inverting the FLR operator, i.e. applying A = 1 - rho^2*nabla_perp^2. struct gkyl_array *flr_rhoSq; // rho^2 weight (times J*g^ij) in the perpendicular Laplacian of A. + struct gkyl_array *flr_kSq; // k^2=1 weight in the FLR operator. + struct gkyl_fem_poisson_perp *flr_op; // Apply the operator 1 - rho_i^2*nabla_perp^2 to invert the FLR operator. struct gkyl_array_integrate *calc_em_energy; // Operator computing EM energy. double *em_energy_red, *em_energy_red_global; // memory for use in GPU reduction of EM energy diff --git a/gyrokinetic/unit/ctest_fem_poisson_perp_lhs_apply.c b/gyrokinetic/unit/ctest_fem_poisson_perp_lhs_apply.c new file mode 100644 index 0000000000..ff63a5b8c8 --- /dev/null +++ b/gyrokinetic/unit/ctest_fem_poisson_perp_lhs_apply.c @@ -0,0 +1,204 @@ +// Test gkyl_fem_poisson_perp_lhs_apply, which applies the LHS operator +// M^{-1}*(M+K) = weak (1 - rho^2*Lap_perp) using the assembled FEM matrices. +// +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +static struct gkyl_array* +mkarr(long nc, long size) +{ + return gkyl_array_new(GKYL_DOUBLE, nc, size); +} + +static struct gkyl_array* +mkarr_dev(bool on_gpu, long nc, long size) +{ +#ifdef GKYL_HAVE_CUDA + if (on_gpu) return gkyl_array_cu_dev_new(GKYL_DOUBLE, nc, size); +#endif + return gkyl_array_new(GKYL_DOUBLE, nc, size); +} + +static double +calc_l2(struct gkyl_rect_grid grid, struct gkyl_range range, struct gkyl_range range_ext, + struct gkyl_basis basis, struct gkyl_array *f1, struct gkyl_array *f2) +{ + struct gkyl_array *diff = mkarr(basis.num_basis, range_ext.volume); + struct gkyl_range_iter iter; + gkyl_range_iter_init(&iter, &range); + while (gkyl_range_iter_next(&iter)) { + long lidx = gkyl_range_idx(&range, iter.idx); + const double *f1p = gkyl_array_cfetch(f1, lidx); + const double *f2p = f2? gkyl_array_cfetch(f2, lidx) : NULL; + double *diffp = gkyl_array_fetch(diff, lidx); + for (int i=0; i 2) { + gkyl_array_shiftc(epsilon_ho, rho*rho*dg0norm, 2*basis.num_basis); + } + struct gkyl_array *kSq_ho = mkarr(basis.num_basis, local_ext.volume); + gkyl_array_shiftc(kSq_ho, -dg0norm, 0); + + struct gkyl_array *src_ho = mkarr(basis.num_basis, local_ext.volume); + gkyl_proj_on_basis *proj = gkyl_proj_on_basis_new(&grid, &basis, poly_order+1, 1, + dim==2? evalFunc_2x : evalFunc_3x, NULL); + gkyl_proj_on_basis_advance(proj, 0.0, &local, src_ho); + gkyl_proj_on_basis_release(proj); + + // Working arrays (on device if use_gpu), populated from the host copies. + struct gkyl_array *epsilon = mkarr_dev(use_gpu, neps*basis.num_basis, local_ext.volume); + struct gkyl_array *kSq = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + struct gkyl_array *src = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + struct gkyl_array *w = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + struct gkyl_array *g = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + struct gkyl_array *w2 = mkarr_dev(use_gpu, basis.num_basis, local_ext.volume); + gkyl_array_copy(epsilon, epsilon_ho); + gkyl_array_copy(kSq, kSq_ho); + gkyl_array_copy(src, src_ho); + + gkyl_fem_poisson_perp *poisson = gkyl_fem_poisson_perp_new(&local, &grid, basis, &bcs, NULL, + epsilon, kSq, use_gpu); + + // Field w = A^{-1}*src, continuous in the perpendicular direction(s). + gkyl_fem_poisson_perp_set_rhs(poisson, src); + gkyl_fem_poisson_perp_solve(poisson, w); + + // Apply the LHS operator, g = (1 - rho^2*Lap_perp) w. + gkyl_fem_poisson_perp_lhs_apply(poisson, w, g); + + // Exact identity: solving with the same operator inverts the apply. + gkyl_fem_poisson_perp_set_rhs(poisson, g); + gkyl_fem_poisson_perp_solve(poisson, w2); + + // Copy results to host for the L2 checks. + struct gkyl_array *w_ho = mkarr(basis.num_basis, local_ext.volume); + struct gkyl_array *g_ho = mkarr(basis.num_basis, local_ext.volume); + struct gkyl_array *w2_ho = mkarr(basis.num_basis, local_ext.volume); + gkyl_array_copy(w_ho, w); + gkyl_array_copy(g_ho, g); + gkyl_array_copy(w2_ho, w2); + + double err = calc_l2(grid, local, local_ext, basis, w2_ho, w_ho); + double norm_w = calc_l2(grid, local, local_ext, basis, w_ho, NULL); + TEST_CHECK(err < 1.0e-10*norm_w); + TEST_MSG("round-trip L2 error = %g (|w| = %g)", err, norm_w); + + // Applying (1 - rho^2*Lap_perp) amplifies the mode, so |g| > |w|. + double norm_g = calc_l2(grid, local, local_ext, basis, g_ho, NULL); + TEST_CHECK(norm_g > norm_w); + TEST_MSG("|g| = %g, |w| = %g", norm_g, norm_w); + + gkyl_array_release(epsilon_ho); + gkyl_array_release(kSq_ho); + gkyl_array_release(src_ho); + gkyl_array_release(epsilon); + gkyl_array_release(kSq); + gkyl_array_release(src); + gkyl_array_release(w); + gkyl_array_release(g); + gkyl_array_release(w2); + gkyl_array_release(w_ho); + gkyl_array_release(g_ho); + gkyl_array_release(w2_ho); + gkyl_fem_poisson_perp_release(poisson); +} + +void test_2x_p1(void) { + int cells[] = {24, 8}; + test_perp_lhs_apply(2, cells, 0.3, GKYL_POISSON_DIRICHLET, false); + test_perp_lhs_apply(2, cells, 0.3, GKYL_POISSON_PERIODIC, false); +} + +void test_3x_p1(void) { + int cells[] = {16, 16, 8}; + test_perp_lhs_apply(3, cells, 0.3, GKYL_POISSON_DIRICHLET, false); + test_perp_lhs_apply(3, cells, 0.3, GKYL_POISSON_PERIODIC, false); +} + +#ifdef GKYL_HAVE_CUDA + +void gpu_test_2x_p1(void) { + int cells[] = {24, 8}; + test_perp_lhs_apply(2, cells, 0.3, GKYL_POISSON_DIRICHLET, true); + test_perp_lhs_apply(2, cells, 0.3, GKYL_POISSON_PERIODIC, true); +} + +void gpu_test_3x_p1(void) { + int cells[] = {16, 16, 8}; + test_perp_lhs_apply(3, cells, 0.3, GKYL_POISSON_DIRICHLET, true); + test_perp_lhs_apply(3, cells, 0.3, GKYL_POISSON_PERIODIC, true); +} + +#endif + +TEST_LIST = { + { "test_2x_p1", test_2x_p1 }, + { "test_3x_p1", test_3x_p1 }, +#ifdef GKYL_HAVE_CUDA + { "gpu_test_2x_p1", gpu_test_2x_p1 }, + { "gpu_test_3x_p1", gpu_test_3x_p1 }, +#endif + { NULL, NULL }, +}; diff --git a/gyrokinetic/zero/fem_poisson_perp.c b/gyrokinetic/zero/fem_poisson_perp.c index ae2e41d20b..3c52312110 100644 --- a/gyrokinetic/zero/fem_poisson_perp.c +++ b/gyrokinetic/zero/fem_poisson_perp.c @@ -76,6 +76,8 @@ gkyl_fem_poisson_perp_new(const struct gkyl_range *solve_range, const struct gky struct gkyl_fem_poisson_perp *up = gkyl_malloc(sizeof(struct gkyl_fem_poisson_perp)); + up->has_lhs_apply = false; + up->bcs = *bcs; up->solve_range = solve_range; up->ndim = grid->ndim; up->ndim_perp = up->ndim-1; @@ -652,8 +654,90 @@ gkyl_fem_poisson_perp_update_lhs(gkyl_fem_poisson_perp *up, struct gkyl_array *e gkyl_superlu_amat_update_from_triples(up->prob, up->tri); } +static void +fem_poisson_perp_lhs_apply_init(gkyl_fem_poisson_perp *up) +{ + // Build a mass-matrix (M) solver: same operator with epsilon=0 and a constant + // kSq=-1 (so the LHS stencil reduces to the mass matrix). + struct gkyl_array *eps_zero, *kSq_mass; +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) { + eps_zero = gkyl_array_cu_dev_new(GKYL_DOUBLE, up->epsilon->ncomp, up->epsilon->size); + kSq_mass = gkyl_array_cu_dev_new(GKYL_DOUBLE, up->num_basis, up->epsilon->size); + } else { + eps_zero = gkyl_array_new(GKYL_DOUBLE, up->epsilon->ncomp, up->epsilon->size); + kSq_mass = gkyl_array_new(GKYL_DOUBLE, up->num_basis, up->epsilon->size); + } +#else + eps_zero = gkyl_array_new(GKYL_DOUBLE, up->epsilon->ncomp, up->epsilon->size); + kSq_mass = gkyl_array_new(GKYL_DOUBLE, up->num_basis, up->epsilon->size); +#endif + gkyl_array_clear(eps_zero, 0.0); + gkyl_array_clear(kSq_mass, 0.0); + gkyl_array_shiftc(kSq_mass, -pow(sqrt(2.0), up->ndim), 0); + + up->mass = gkyl_fem_poisson_perp_new(up->solve_range, &up->grid, up->basis, &up->bcs, NULL, + eps_zero, kSq_mass, up->use_gpu); + + long numnodes_tot = up->numnodes_global*up->par_range.volume; + up->lhs_dual = gkyl_malloc(sizeof(double[numnodes_tot])); +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) + up->lhs_dual_cu = gkyl_cu_malloc(sizeof(double[numnodes_tot])); +#endif + + gkyl_array_release(eps_zero); + gkyl_array_release(kSq_mass); + up->has_lhs_apply = true; +} + +void +gkyl_fem_poisson_perp_lhs_apply(gkyl_fem_poisson_perp *up, struct gkyl_array *xin, struct gkyl_array *xout) +{ + assert(up->ishelmholtz); + assert(up->num_bias_line == 0); // Biased mass solve not implemented. + if (!up->has_lhs_apply) fem_poisson_perp_lhs_apply_init(up); + +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) { + // Recover the global nodal vector x_nodal = M^{-1}*(M_src*xin). + gkyl_fem_poisson_perp_set_rhs(up->mass, xin); + gkyl_culinsolver_solve(up->mass->prob_cu); + gkyl_culinsolver_sync(up->mass->prob_cu); + double *x_nodal = gkyl_culinsolver_get_sol_ptr(up->mass->prob_cu, 0); + + // Dual d = (M+K)*x_nodal, then solve M*z = d (-> modal xout). + gkyl_culinsolver_mat_vec(up->prob_cu, x_nodal, up->lhs_dual_cu); + double *rhs_cu = gkyl_culinsolver_get_rhs_ptr(up->mass->prob_cu, 0); + gkyl_cu_memcpy(rhs_cu, up->lhs_dual_cu, + sizeof(double)*up->numnodes_global*up->par_range.volume, GKYL_CU_MEMCPY_D2D); + gkyl_fem_poisson_perp_solve(up->mass, xout); + return; + } +#endif + + // Recover the global nodal vector x_nodal = M^{-1}*(M_src*xin). + gkyl_fem_poisson_perp_set_rhs(up->mass, xin); + gkyl_superlu_solve(up->mass->prob); + double *x_nodal = gkyl_superlu_get_rhs_ptr(up->mass->prob, 0); + + // Dual d = (M+K)*x_nodal, then solve M*z = d (-> modal xout). + gkyl_superlu_mat_vec(up->prob, x_nodal, up->lhs_dual); + gkyl_superlu_brhs_from_array(up->mass->prob, up->lhs_dual); + gkyl_fem_poisson_perp_solve(up->mass, xout); +} + void gkyl_fem_poisson_perp_release(struct gkyl_fem_poisson_perp *up) { + if (up->has_lhs_apply) { + gkyl_fem_poisson_perp_release(up->mass); + gkyl_free(up->lhs_dual); +#ifdef GKYL_HAVE_CUDA + if (up->use_gpu) + gkyl_cu_free(up->lhs_dual_cu); +#endif + } + if (up->isdomperiodic) { gkyl_array_release(up->rhs_cellavg); gkyl_free(up->rhs_avg); diff --git a/gyrokinetic/zero/gkyl_fem_poisson_perp.h b/gyrokinetic/zero/gkyl_fem_poisson_perp.h index d68db9e9f8..bf5508dc3e 100644 --- a/gyrokinetic/zero/gkyl_fem_poisson_perp.h +++ b/gyrokinetic/zero/gkyl_fem_poisson_perp.h @@ -65,6 +65,19 @@ void gkyl_fem_poisson_perp_solve(gkyl_fem_poisson_perp* up, struct gkyl_array *p */ void gkyl_fem_poisson_perp_update_lhs(gkyl_fem_poisson_perp* up, struct gkyl_array *epsilon, struct gkyl_array *kSq); +/** + * Apply the LHS operator, i.e. compute xout = M^{-1}*(M+K)*xin where M is the + * mass matrix and M+K is the (perpendicular Helmholtz) LHS matrix. With + * epsilon=rho^2 and kSq=-1 this returns the weak image of + * (1 - rho^2*Lap_perp)*xin. Assumes xin is continuous in the perpendicular + * directions and requires a Helmholtz solver (kSq!=NULL). + * + * @param up FEM poisson updater to run. + * @param xin DG field to apply the operator to. + * @param xout DG field holding the result. + */ +void gkyl_fem_poisson_perp_lhs_apply(gkyl_fem_poisson_perp* up, struct gkyl_array *xin, struct gkyl_array *xout); + /** * Delete updater. * diff --git a/gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h b/gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h index 3ed0c771d3..051fe2d476 100644 --- a/gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h +++ b/gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h @@ -621,6 +621,13 @@ struct gkyl_fem_poisson_perp { int num_bias_line; // Number of biased lines. struct gkyl_poisson_bias_line *bias_lines; // Biased lines. bias_src_func_t bias_line_src; // Function to enforce biasing in RHS source. + + // Objects to apply the LHS operator (M+K), used to apply 1-rho^2*Lap_perp. + bool has_lhs_apply; // Whether the LHS-apply objects are allocated. + struct gkyl_poisson_bc bcs; // Saved BCs (used to build the mass solver). + struct gkyl_fem_poisson_perp *mass; // Mass-matrix (M) solver for the LHS apply. + double *lhs_dual; // Host global dual (M+K)*x vector. + double *lhs_dual_cu; // Device global dual (M+K)*x vector. }; void From 41216a8313ee045d9631f0c1d474bb609b8cc5d7 Mon Sep 17 00:00:00 2001 From: Antoinehoff Date: Mon, 6 Jul 2026 07:07:45 -0700 Subject: [PATCH 04/12] simplify some comments --- gyrokinetic/apps/gk_field_flr.c | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/gyrokinetic/apps/gk_field_flr.c b/gyrokinetic/apps/gk_field_flr.c index 2cbc164db5..26f61a399f 100644 --- a/gyrokinetic/apps/gk_field_flr.c +++ b/gyrokinetic/apps/gk_field_flr.c @@ -22,12 +22,8 @@ gk_field_flr_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f) rhoSq_ref = pow(f->info.flr.gyroradius, 2.0); } else { - // GKYL_GK_FLR_PADE_CONST_SUM, or field-level FLR options not set (default): // polarization-weighted average of the species gyroradii, // rho^2 = sum_s eps_s0*rho_s0^2 / sum_s eps_s0, eps_s0 = n_s0*m_s/B^2, - // which suppresses the electron contribution by m_e/m_i and reduces to - // rho_i0^2 for a single ion species. Species without FLR enabled - // contribute rho_s0^2 = 0. double polarization_bmag = f->info.polarization_bmag ? f->info.polarization_bmag : app->bmag_ref; double eps_sum = 0.0; for (int i = 0; i < app->num_species; ++i) { @@ -40,8 +36,7 @@ gk_field_flr_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f) } rhoSq_ref /= eps_sum; } - // Weight in the perpendicular Laplacian of A, staged for the forthcoming - // Laplacian-apply infrastructure. + // Weight in the perpendicular Laplacian. f->flr_rhoSq = mkarr(app->use_gpu, (2 * (app->cdim - 1) - 1) * app->basis.num_basis, app->local_ext.volume); gkyl_array_set_offset(f->flr_rhoSq, rhoSq_ref, app->gk_geom->geo_int.gxxj, 0 * app->basis.num_basis); if (app->cdim > 2) { From 73a17bd6fa8937775c825c581f226c992f81534d Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Mon, 6 Jul 2026 10:21:07 -0400 Subject: [PATCH 05/12] add a regression test with FLR effects. --- gyrokinetic/creg/rt_gk_flr_d3d_iwl_2x2v_p1.c | 836 +++++++++++++++++++ 1 file changed, 836 insertions(+) create mode 100644 gyrokinetic/creg/rt_gk_flr_d3d_iwl_2x2v_p1.c diff --git a/gyrokinetic/creg/rt_gk_flr_d3d_iwl_2x2v_p1.c b/gyrokinetic/creg/rt_gk_flr_d3d_iwl_2x2v_p1.c new file mode 100644 index 0000000000..399419014d --- /dev/null +++ b/gyrokinetic/creg/rt_gk_flr_d3d_iwl_2x2v_p1.c @@ -0,0 +1,836 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +// Define the context of the simulation. This is basically all the globals +struct gk_app_ctx { + int cdim, vdim; // Dimensionality. + + // Geometry and magnetic field. + double a_shift; // Parameter in Shafranov shift. + double Z_axis; // Magnetic axis height [m]. + double R_axis; // Magnetic axis major radius [m]. + + double R0; // Major radius of the simulation box [m]. + double a_mid; // Minor radius at outboard midplane [m]. + double r0; // Minor radius of the simulation box [m]. + double B0; // Magnetic field magnitude in the simulation box [T]. + double kappa; // Elongation (=1 for no elongation). + double delta; // Triangularity (=0 for no triangularity). + double q0; // Magnetic safety factor in the center of domain. + double Cy; // Prefactor in binormal coordinate. + + double x_LCFS; // Radial location of the last closed flux surface. + + // Plasma parameters. + double me; double qe; + double mi; double qi; + double n0; double Te0; double Ti0; + double rho_i; + + // Collisions. + double nu_frac; + + // Source parameters. + double n_srcOMP; // Amplitude of the OMP source + double x_srcOMP; // Radial location of the OMP source. + double Te_srcOMP; // Te for the OMP source. + double Ti_srcOMP; // Ti for the OMP source. + double sigma_srcOMP; // Radial spread of the OMP source. + double n_srcGB; // Amplitude of the grad-B source + double x_srcGB; // Radial location of the grad-B source. + double sigma_srcGB; // Radial spread of the grad-B source. + double bfac_srcGB; // Field aligned spread of the grad-B source. + double Te_srcGB; // Te for the grad-B source. + double Ti_srcGB; // Ti for the grad-B source. + double floor_src; // Source floor. + + // Grid parameters. + double Lx; // Domain size in radial direction. + double Lz; // Domain size along magnetic field. + double x_min; double x_max; + double z_min; double z_max; + int Nx; + int Nz; + int Nvpar; + int Nmu; + int cells[GKYL_MAX_DIM]; // Number of cells in all directions. + int poly_order; + double vpar_max_elc; double mu_max_elc; + double vpar_max_ion; double mu_max_ion; + + double t_end; // End time. + int num_frames; // Number of output frames. + double write_phase_freq; // Frequency of writing phase-space diagnostics (as a fraction of num_frames). + int int_diag_calc_num; // Number of integrated diagnostics computations (=INT_MAX for every step). + double dt_failure_tol; // Minimum allowable fraction of initial time-step. + int num_failures_max; // Maximum allowable number of consecutive small time-steps. +}; + +double r_x(double x, double a_mid) +{ + return x+a_mid-0.1; +} + +double qprofile(double R) +{ + // Magnetic safety factor as a function of minor radius r. + double a[] = {384.2108662176567, -2417.1050263480997, 5079.5143390131225, -3563.038325842346}; + return a[0]*pow(R,3) + a[1]*pow(R,2) + a[2]*(R) + a[3]; +} + +double R_rtheta(double r, double theta, void *ctx) +{ + // Major radius as a function of minor radius r and poloidal angle theta. + struct gk_app_ctx *app = ctx; + double a_shift = app->a_shift; + double R_axis = app->R_axis; + double delta = app->delta; + return R_axis - a_shift*r*r/(2.*R_axis) + r*cos(theta + asin(delta)*sin(theta)); +} + +double Z_rtheta(double r, double theta, void *ctx) +{ + // Z (height) as a function of minor radius r and poloidal angle theta. + struct gk_app_ctx *app = ctx; + double Z_axis = app->Z_axis; + double kappa = app->kappa; + return Z_axis + kappa*r*sin(theta); +} + +// Partial derivatives of R(r,theta) and Z(r,theta) +double dRdr(double r, double theta, void *ctx) +{ + struct gk_app_ctx *app = ctx; + double a_shift = app->a_shift; + double R_axis = app->R_axis; + double delta = app->delta; + return - a_shift*r/(R_axis) + cos(theta + asin(delta)*sin(theta)); +} +double dRdtheta(double r, double theta, void *ctx) +{ + struct gk_app_ctx *app = ctx; + double delta = app->delta; + return -r*sin(theta + asin(delta)*sin(theta))*(1.+asin(delta)*cos(theta)); +} +double dZdr(double r, double theta, void *ctx) +{ + struct gk_app_ctx *app = ctx; + double kappa = app->kappa; + return kappa*sin(theta); +} +double dZdtheta(double r, double theta, void *ctx) +{ + struct gk_app_ctx *app = ctx; + double kappa = app->kappa; + return kappa*r*cos(theta); +} + +double Jr(double r, double theta, void *ctx) +{ + return R_rtheta(r,theta,ctx)*( dRdr(r,theta,ctx) * dZdtheta(r,theta,ctx) + -dRdtheta(r,theta,ctx) * dZdr(r,theta,ctx) ); +} + +struct integrand_ctx { + struct gk_app_ctx *app_ctx; + double r; +}; +double integrand(double t, void *int_ctx) +{ + struct integrand_ctx *inctx = int_ctx; + double r = inctx->r; + struct gk_app_ctx *app = inctx->app_ctx; + return Jr(r,t,app) / pow(R_rtheta(r,t,app),2); +} + +double Bphi(double R, void *ctx) +{ + // Toroidal magnetic field. + struct gk_app_ctx *app = ctx; + double B0 = app->B0; + double R0 = app->R0; + return B0*R0/R; +} + +double dPsidr(double r, double theta, void *ctx) +{ + struct gk_app_ctx *app = ctx; + struct integrand_ctx tmp_ctx = {.app_ctx = app, .r = r}; + struct gkyl_qr_res integral; + integral = gkyl_dbl_exp(integrand, &tmp_ctx, 0., 2.*M_PI, 7, 1e-10); + + double R = R_rtheta(r,theta,ctx); + double Bt = Bphi(R,ctx); + double R_omp = R_rtheta(r,0.0,ctx); + return ( R*Bt/(2.*M_PI*qprofile(R_omp)) )*integral.res; +} + +double alpha(double r, double theta, double phi, void *ctx) +{ + double twrap = theta; + while (twrap < -M_PI) twrap = twrap+2.*M_PI; + while (M_PI < twrap) twrap = twrap-2.*M_PI; + + struct gk_app_ctx *app = ctx; + struct integrand_ctx tmp_ctx = {.app_ctx = app, .r = r}; + struct gkyl_qr_res integral; + if (0. < twrap) { + integral = gkyl_dbl_exp(integrand, &tmp_ctx, 0., twrap, 7, 1e-10); + } else { + integral = gkyl_dbl_exp(integrand, &tmp_ctx, twrap, 0., 7, 1e-10); + integral.res = -integral.res; + } + + double R = R_rtheta(r,theta,ctx); + double Bt = Bphi(R,ctx); + + return phi - R*Bt*integral.res/dPsidr(r,theta,ctx); +} + +double gradr(double r, double theta, void *ctx) +{ + return (R_rtheta(r,theta,ctx)/Jr(r,theta,ctx))*sqrt(pow(dRdtheta(r,theta,ctx),2) + pow(dZdtheta(r,theta,ctx),2)); +} + +// Common source profiles. +void density_srcOMP(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double n_srcOMP = app->n_srcOMP; + double x_srcOMP = app->x_srcOMP; + double sigma_srcOMP = app->sigma_srcOMP; + double floor_src = app->floor_src; + + fout[0] = n_srcOMP*(exp(-(pow(x-x_srcOMP,2))/(2.*pow(sigma_srcOMP,2)))+floor_src); +} +void zero_func(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + fout[0] = 0.0; +} + +// Electron source profiles. +void density_elc_srcGB(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double n_srcGB = app->n_srcGB; + double x_srcGB = app->x_srcGB; + double sigma_srcGB = app->sigma_srcGB; + double bfac_srcGB = app->bfac_srcGB; + double floor_src = app->floor_src; + + fout[0] = n_srcGB*exp(-pow(x-x_srcGB,2)/(2.*pow(sigma_srcGB,2))) + *GKYL_MAX2(sin(z)*exp(-pow(fabs(z),1.5)/(2*pow(bfac_srcGB,2))),0.); +} +void temp_elc_srcOMP(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double x_srcOMP = app->x_srcOMP; + double sigma_srcOMP = app->sigma_srcOMP; + double Te_srcOMP = app->Te_srcOMP; + + if (x < x_srcOMP + 3*sigma_srcOMP) { + fout[0] = Te_srcOMP; + } else { + fout[0] = Te_srcOMP*3./8.; + } +} +void temp_elc_srcGB(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double Te_srcGB = app->Te_srcGB; + + fout[0] = Te_srcGB; +} + +// Ion source profiles. +void density_ion_srcGB(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double n_srcGB = app->n_srcGB; + double x_srcGB = app->x_srcGB; + double sigma_srcGB = app->sigma_srcGB; + double bfac_srcGB = app->bfac_srcGB; + double floor_src = app->floor_src; + + fout[0] = -n_srcGB*exp(-pow(x-x_srcGB,2)/(2.*pow(sigma_srcGB,2))) + *GKYL_MAX2(sin(z)*exp(-pow(fabs(z),1.5)/(2*pow(bfac_srcGB,2))),0.); +} +void temp_ion_srcOMP(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double x_srcOMP = app->x_srcOMP; + double sigma_srcOMP = app->sigma_srcOMP; + double Ti_srcOMP = app->Ti_srcOMP; + + if (x < x_srcOMP + 3*sigma_srcOMP) { + fout[0] = Ti_srcOMP; + } else { + fout[0] = Ti_srcOMP*3./8.; + } +} +void temp_ion_srcGB(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double Ti_srcGB = app->Ti_srcGB; + + fout[0] = Ti_srcGB; +} + +// Ion initial conditions +void density_init(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double n0 = app->n0; + + fout[0] = n0*1.1*(0.5*(1.+tanh(2.*(2.-25.*x)))+0.01); +} +void temp_ion(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double Ti0 = app->Ti0; + + fout[0] = Ti0*((1./3.)*(2.+tanh(2.*(2.-25.*x)))+0.01); +} + +// Electron initial conditions +void temp_elc(double t, const double * GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void *ctx) +{ + double x = xn[0], z = xn[1]; + + struct gk_app_ctx *app = ctx; + double Te0 = app->Te0; + + fout[0] = Te0*((1./3.)*(2.+tanh(2.*(2.-25.*x)))+0.01); +} + +void +diffusion_D_func(double t, const double* GKYL_RESTRICT xn, double* GKYL_RESTRICT fout, void* ctx) +{ + struct sheath_ctx *app = ctx; + + fout[0] = 0.5; // Diffusivity [m^2/s]. +} + +// Geometry evaluation functions for the gk app +void mapc2p(double t, const double *xc, double* GKYL_RESTRICT xp, void *ctx) +{ + struct gk_app_ctx *app = ctx; + double Cy = app->Cy; + double a_mid = app->a_mid; + + double x = xc[0], y = xc[1], z = xc[2]; + + double r = r_x(x,a_mid); + + // Map to cylindrical (R, Z, phi) coordinates. + double R = R_rtheta(r, z, ctx); + double Z = Z_rtheta(r, z, ctx); + double phi = y/Cy + alpha(r, z, 0, ctx); + // Map to Cartesian (X, Y, Z) coordinates. + double X = R*cos(phi); + double Y = R*sin(phi); + + xp[0] = X; xp[1] = Y; xp[2] = Z; +} + +void bfield_func(double t, const double *xc, double* GKYL_RESTRICT fout, void *ctx) +{ + struct gk_app_ctx *app = ctx; + double Cy = app->Cy; + double a_mid = app->a_mid; + + double x = xc[0], y = xc[1], z = xc[2]; + double r = r_x(x,a_mid); + double Bt = Bphi(R_rtheta(r,z,ctx),ctx); + double Bp = dPsidr(r,z,ctx)/R_rtheta(r,z,ctx)*gradr(r,z,ctx); + + double drdtheta = dRdtheta(r,z,ctx); + double dzdtheta = dZdtheta(r,z,ctx); + double den = sqrt(pow(drdtheta,2) + pow(dzdtheta,2)); + double B_r = Bp*drdtheta/den; + double B_z = Bp*dzdtheta/den; + double phi = y/Cy + alpha(r, z, 0, ctx); + double R = R_rtheta(r, z, ctx); + + // xc are computational coords. + // Set Cartesian components of magnetic field. + fout[0] = B_r * cos(phi) + Bt * sin(phi); + fout[1] = B_r * sin(phi) - Bt * cos(phi); + fout[2] = B_z; +} + +struct gk_app_ctx +create_ctx(void) +{ + int cdim = 2, vdim = 2; // Dimensionality. + + // Universal constant parameters. + double eps0 = GKYL_EPSILON0, eV = GKYL_ELEMENTARY_CHARGE; + double mp = GKYL_PROTON_MASS, me = GKYL_ELECTRON_MASS; + double qi = eV; // ion charge + double qe = -eV; // electron charge + + // Geometry and magnetic field. + double a_shift = 0.6; // Parameter in Shafranov shift. + double Z_axis = 0.00232616113; // Magnetic axis height [m]. + double R_axis = 1.72068012; // Magnetic axis major radius [m]. + double B_axis = 2.0; // Magnetic field at the magnetic axis [T]. + double R_LCFSmid = 2.2801477223421736; // Major radius of the LCFS at the outboard midplane [m]. + double Rmid_min = R_LCFSmid - 5*0.15/8; // Minimum midplane major radius of simulation box [m]. + double Rmid_max = R_LCFSmid + 3*0.15/8; // Maximum midplane major radius of simulation box [m]. + double R0 = 0.5*(Rmid_min+Rmid_max); // Major radius of the simulation box [m]. + + double a_mid = R_LCFSmid-R_axis; // Minor radius at outboard midplane [m]. + // Redefine a_mid with Shafranov shift, to ensure LCFS radial location. + a_mid = R_axis/a_shift - sqrt(R_axis*(R_axis - 2*a_shift*R_LCFSmid + 2*a_shift*R_axis))/a_shift; + + double r0 = R0-R_axis; // Minor radius of the simulation box [m]. + double B0 = B_axis*(R_axis/R0); // Magnetic field magnitude in the simulation box [T]. + double kappa = 1.488; // Elongation (=1 for no elongation). + double delta = 0.1; // Triangularity (=0 for no triangularity). + double Lx = Rmid_max-Rmid_min; // Domain size along x. + double Lz = 2.*M_PI-1e-10; // Domain size along magnetic field. + double x_min = 0.; + double x_max = Lx; + double z_min = -Lz/2.; + double z_max = Lz/2.; + + double x_LCFS = R_LCFSmid - Rmid_min; // Radial location of the last closed flux surface. + + double q0 = qprofile(R0); // Magnetic safety factor in the center of domain. + double Cy = 1; // Normalization in binormal coordinate. + + // Plasma parameters. Chosen based on the value of a cubic sline + // between the last TS data inside the LCFS and the probe data in + // in the far SOL, near R=0.475 m. + double AMU = 2.01410177811; + double mi = mp*AMU; // Deuterium ions. + double Te0 = 200*eV; + double Ti0 = 200*eV; + double n0 = 2.0e19; // [1/m^3] + + double nu_frac = 1.0; + + double vte = sqrt(Te0/me), vti = sqrt(Ti0/mi); // Thermal speeds. + + double c_s = sqrt(Te0/mi); + double omega_ci = fabs(qi*B0/mi); + double rho_s = c_s/omega_ci; + double rho_i = vti/omega_ci; + + // Source parameters + double n_srcOMP = 3.e19; + double x_srcOMP = x_min; + double Te_srcOMP = 2*Te0; + double Ti_srcOMP = 2*Ti0; + double sigma_srcOMP = 0.03*Lx; + double n_srcGB = 0.04*1.3553306021427228e+23; + double x_srcGB = x_min; + double sigma_srcGB = 10*rho_s; + double bfac_srcGB = 1.2; + double Te_srcGB = 350*eV; + double Ti_srcGB = 350*eV; + double floor_src = 1e-2; + + // Grid parameters + int Nx = 8; + int Nz = 12; + int Nvpar = 8; + int Nmu = 4; + int poly_order = 1; + + double vpar_max_elc = 4.*vte; + double mu_max_elc = me*pow(4*vte,2)/(2*B0); + double vpar_max_ion = 4.*vti; + double mu_max_ion = mi*pow(4*vti,2)/(2*B0); + + double t_end = 1.e-7; + int num_frames = 1; + double write_phase_freq = 0.2; // Frequency of writing phase-space diagnostics (as a fraction of num_frames). + int int_diag_calc_num = num_frames*100; + double dt_failure_tol = 1.0e-4; // Minimum allowable fraction of initial time-step. + int num_failures_max = 20; // Maximum allowable number of consecutive small time-steps. + + struct gk_app_ctx ctx = { + .cdim = cdim, + .vdim = vdim, + .a_shift = a_shift, + .R_axis = R_axis, + .R0 = R0 , + .a_mid = a_mid , + .r0 = r0 , + .B0 = B0 , + .kappa = kappa , + .delta = delta , + .q0 = q0 , + .Cy = Cy , + .Lx = Lx , + .Lz = Lz , + .x_min = x_min, .x_max = x_max, + .z_min = z_min, .z_max = z_max, + + .x_LCFS = x_LCFS, + + .me = me, .qe = qe, + .mi = mi, .qi = qi, + .n0 = n0, .Te0 = Te0, .Ti0 = Ti0, + .rho_i = rho_i, + + .nu_frac = nu_frac, + + .n_srcOMP = n_srcOMP , + .x_srcOMP = x_srcOMP , + .Te_srcOMP = Te_srcOMP , + .Ti_srcOMP = Ti_srcOMP , + .sigma_srcOMP = sigma_srcOMP, + .n_srcGB = n_srcGB , + .x_srcGB = x_srcGB , + .sigma_srcGB = sigma_srcGB , + .bfac_srcGB = bfac_srcGB , + .Te_srcGB = Te_srcGB , + .Ti_srcGB = Ti_srcGB , + .floor_src = floor_src , + + .Nx = Nx, + .Nz = Nz, + .Nvpar = Nvpar, + .Nmu = Nmu, + .cells = {Nx, Nz, Nvpar, Nmu}, + .poly_order = poly_order, + .vpar_max_elc = vpar_max_elc, .mu_max_elc = mu_max_elc, + .vpar_max_ion = vpar_max_ion, .mu_max_ion = mu_max_ion, + + .t_end = t_end, .num_frames = num_frames, + .write_phase_freq = write_phase_freq, + .int_diag_calc_num = int_diag_calc_num, + .dt_failure_tol = dt_failure_tol, + .num_failures_max = num_failures_max, + }; + return ctx; +} + + +int +main(int argc, char **argv) +{ + struct gkyl_app_args app_args = parse_app_args(argc, argv); + +#ifdef GKYL_HAVE_MPI + if (app_args.use_mpi) MPI_Init(&argc, &argv); +#endif + + if (app_args.trace_mem) { + gkyl_cu_dev_mem_debug_set(true); + gkyl_mem_debug_set(true); + } + + struct gk_app_ctx ctx = create_ctx(); // Context for init functions. + + int cells_x[ctx.cdim], cells_v[ctx.vdim]; + for (int d=0; d Date: Mon, 6 Jul 2026 11:43:53 -0400 Subject: [PATCH 06/12] revert the change of operator for the species FLR effect. We have to use deflated here because we use the value at the boundary as BC. --- gyrokinetic/apps/gk_species.c | 32 ++++++++++++++++++------ gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 2 +- 2 files changed, 26 insertions(+), 8 deletions(-) diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index fe3acd9905..c59a4aedbd 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/gyrokinetic/apps/gk_species.c @@ -33,8 +33,7 @@ gk_species_gyroaverage_enabled(gkyl_gyrokinetic_app *app, struct gk_species *spe // get the gyroaveraged potential used in the Hamiltonian. struct timespec wst = gkyl_wall_clock(); // Gyroaverage input field. - gkyl_fem_poisson_perp_set_rhs(species->flr_op, field_in); - gkyl_fem_poisson_perp_solve(species->flr_op, field_gyroavg); + gkyl_deflated_fem_poisson_advance(species->flr_op, field_in, field_in, field_gyroavg); app->stat.species_gyroavg_tm += gkyl_time_diff_now_sec(wst); } @@ -1576,10 +1575,29 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st gks->flr_kSq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); gkyl_array_shiftc(gks->flr_kSq, -pow(sqrt(2.0),app->cdim), 0); // Sets kSq=-1. - // Use the same solver and BCs as the field solve so the gyroaverage and - // the FLR inversion of phi use the same discretization. - gks->flr_op = gkyl_fem_poisson_perp_new(&app->local, &app->grid, app->basis, - &app->field->poisson_bcs, 0, gks->flr_rhoSqD2, gks->flr_kSq, app->use_gpu); + // The gyroaverage uses varying-Dirichlet BCs (input field as its own + // boundary value) so no FLR correction is applied at the wall. + struct gkyl_poisson_bc flr_bc; + for (int d=0; dcdim-1; d++) { + struct gkyl_gyrokinetic_bc *bc_lo = gk_fetch_bc_with_dir_edge(app->field->info.poisson_bcs, 2*app->cdim, d, GKYL_LOWER_EDGE); + if (bc_lo != 0) { + if (bc_lo->type == GKYL_BC_GK_FIELD_PERIODIC) + flr_bc.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); + else + flr_bc.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_DIRICHLET_VARYING); + } + + struct gkyl_gyrokinetic_bc *bc_up = gk_fetch_bc_with_dir_edge(app->field->info.poisson_bcs, 2*app->cdim, d, GKYL_UPPER_EDGE); + if (bc_up != 0) { + if (bc_up->type == GKYL_BC_GK_FIELD_PERIODIC) + flr_bc.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); + else + flr_bc.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_DIRICHLET_VARYING); + } + } + // Deflated Poisson solve is performed on range assuming decomposition is *only* in z. + gks->flr_op = gkyl_deflated_fem_poisson_new(app->grid, app->basis_on_dev, app->basis, + app->local, app->local, gks->flr_rhoSqD2, gks->flr_kSq, flr_bc, NULL, app->use_gpu); } else { gks->gyroaverage = gk_species_gyroaverage_disabled; @@ -2045,7 +2063,7 @@ gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *gks if (gks->info.flr.type) { gkyl_array_release(gks->flr_rhoSqD2); gkyl_array_release(gks->flr_kSq); - gkyl_fem_poisson_perp_release(gks->flr_op); + gkyl_deflated_fem_poisson_release(gks->flr_op); } gks->release_func(app, gks); diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index a300d749a4..7c0eb15afd 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -1148,7 +1148,7 @@ struct gk_species { struct gkyl_array *m0_gyroavg; // Gyroaveraged particle density. struct gkyl_array *flr_rhoSqD2; // Laplacian weight in FLR operator. struct gkyl_array *flr_kSq; // Field multiplying phi in FLR operator. - struct gkyl_fem_poisson_perp *flr_op; // Screened Poisson solver used to gyroaverage fields. + struct gkyl_deflated_fem_poisson *flr_op; // Screened Poisson solver used to gyroaverage fields. // Pointer to function that performs the gyroaverage. void (*gyroaverage)(gkyl_gyrokinetic_app *app, struct gk_species *species, struct gkyl_array *field_in, struct gkyl_array *field_gyroavg); From c32b8a7018e037a9722e7a2f6a6a75db5e2434fe Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Mon, 6 Jul 2026 11:59:28 -0400 Subject: [PATCH 07/12] add a copy of the bias list to the host so that we can pass it to the apply LHS builder and build the exact same matrix. I am not sure on how stable and sound this is but otherwise, the FLR operator will break the bias plane BC. In general, I am not sure that the clopen setup is compatible with FLR since this is varying BCs... --- .../unit/ctest_fem_poisson_perp_lhs_apply.c | 41 ++++++++++++++----- gyrokinetic/zero/fem_poisson_perp.c | 20 ++++++++- gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h | 2 + moments/zero/fem_poisson.c | 19 ++++++++- moments/zero/gkyl_fem_poisson_priv.h | 2 + 5 files changed, 71 insertions(+), 13 deletions(-) diff --git a/gyrokinetic/unit/ctest_fem_poisson_perp_lhs_apply.c b/gyrokinetic/unit/ctest_fem_poisson_perp_lhs_apply.c index ff63a5b8c8..afb85882c2 100644 --- a/gyrokinetic/unit/ctest_fem_poisson_perp_lhs_apply.c +++ b/gyrokinetic/unit/ctest_fem_poisson_perp_lhs_apply.c @@ -68,7 +68,8 @@ static void evalFunc_3x(double t, const double *xn, double *fout, void *ctx) } static void -test_perp_lhs_apply(int dim, int cells[], double rho, enum gkyl_poisson_bc_type bc_type, bool use_gpu) +test_perp_lhs_apply(int dim, int cells[], double rho, enum gkyl_poisson_bc_type bc_type, + struct gkyl_poisson_bias_line_list *bias_list, bool use_gpu) { double lower[GKYL_MAX_CDIM], upper[GKYL_MAX_CDIM]; for (int d=0; dhas_lhs_apply = false; up->bcs = *bcs; + // Keep a host-side copy of the bias line list (used to build the mass + // solver in the LHS apply). + up->bias_list_ho.num_bias_line = 0; + up->bias_list_ho.bl = 0; + if (bias_lines) { + if (bias_lines->num_bias_line > 0) { + up->bias_list_ho.num_bias_line = bias_lines->num_bias_line; + size_t bl_sz = bias_lines->num_bias_line * sizeof(struct gkyl_poisson_bias_line); + up->bias_list_ho.bl = gkyl_malloc(bl_sz); + memcpy(up->bias_list_ho.bl, bias_lines->bl, bl_sz); + } + } up->solve_range = solve_range; up->ndim = grid->ndim; up->ndim_perp = up->ndim-1; @@ -676,7 +688,9 @@ fem_poisson_perp_lhs_apply_init(gkyl_fem_poisson_perp *up) gkyl_array_clear(kSq_mass, 0.0); gkyl_array_shiftc(kSq_mass, -pow(sqrt(2.0), up->ndim), 0); - up->mass = gkyl_fem_poisson_perp_new(up->solve_range, &up->grid, up->basis, &up->bcs, NULL, + // Bias the mass solver like the LHS so biased values pass through the apply. + struct gkyl_poisson_bias_line_list *bias_list = up->bias_list_ho.num_bias_line > 0? &up->bias_list_ho : NULL; + up->mass = gkyl_fem_poisson_perp_new(up->solve_range, &up->grid, up->basis, &up->bcs, bias_list, eps_zero, kSq_mass, up->use_gpu); long numnodes_tot = up->numnodes_global*up->par_range.volume; @@ -695,7 +709,6 @@ void gkyl_fem_poisson_perp_lhs_apply(gkyl_fem_poisson_perp *up, struct gkyl_array *xin, struct gkyl_array *xout) { assert(up->ishelmholtz); - assert(up->num_bias_line == 0); // Biased mass solve not implemented. if (!up->has_lhs_apply) fem_poisson_perp_lhs_apply_init(up); #ifdef GKYL_HAVE_CUDA @@ -738,6 +751,9 @@ void gkyl_fem_poisson_perp_release(struct gkyl_fem_poisson_perp *up) #endif } + if (up->bias_list_ho.num_bias_line > 0) + gkyl_free(up->bias_list_ho.bl); + if (up->isdomperiodic) { gkyl_array_release(up->rhs_cellavg); gkyl_free(up->rhs_avg); diff --git a/gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h b/gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h index 051fe2d476..4002249231 100644 --- a/gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h +++ b/gyrokinetic/zero/gkyl_fem_poisson_perp_priv.h @@ -625,6 +625,8 @@ struct gkyl_fem_poisson_perp { // Objects to apply the LHS operator (M+K), used to apply 1-rho^2*Lap_perp. bool has_lhs_apply; // Whether the LHS-apply objects are allocated. struct gkyl_poisson_bc bcs; // Saved BCs (used to build the mass solver). + struct gkyl_poisson_bias_line_list bias_list_ho; // Host-side copy of the bias + // line list (used to build the mass solver). struct gkyl_fem_poisson_perp *mass; // Mass-matrix (M) solver for the LHS apply. double *lhs_dual; // Host global dual (M+K)*x vector. double *lhs_dual_cu; // Device global dual (M+K)*x vector. diff --git a/moments/zero/fem_poisson.c b/moments/zero/fem_poisson.c index 01941a72c6..ffe1143c00 100644 --- a/moments/zero/fem_poisson.c +++ b/moments/zero/fem_poisson.c @@ -64,6 +64,18 @@ gkyl_fem_poisson_new(const struct gkyl_range *solve_range, const struct gkyl_rec up->has_lhs_apply = false; up->bcs = *bcs; + // Keep a host-side copy of the bias plane list (used to build the mass + // solver in the LHS apply). + up->bias_list_ho.num_bias_plane = 0; + up->bias_list_ho.bp = 0; + if (bias_planes) { + if (bias_planes->num_bias_plane > 0) { + up->bias_list_ho.num_bias_plane = bias_planes->num_bias_plane; + size_t bp_sz = bias_planes->num_bias_plane * sizeof(struct gkyl_poisson_bias_plane); + up->bias_list_ho.bp = gkyl_malloc(bp_sz); + memcpy(up->bias_list_ho.bp, bias_planes->bp, bp_sz); + } + } up->solve_range = solve_range; up->ndim = grid->ndim; up->grid = *grid; @@ -428,7 +440,9 @@ fem_poisson_lhs_apply_init(gkyl_fem_poisson *up) gkyl_array_clear(kSq_mass, 0.0); gkyl_array_shiftc(kSq_mass, -pow(sqrt(2.0), up->ndim), 0); - up->mass = gkyl_fem_poisson_new(up->solve_range, &up->grid, up->basis, &up->bcs, NULL, + // Bias the mass solver like the LHS so biased values pass through the apply. + struct gkyl_poisson_bias_plane_list *bias_list = up->bias_list_ho.num_bias_plane > 0? &up->bias_list_ho : NULL; + up->mass = gkyl_fem_poisson_new(up->solve_range, &up->grid, up->basis, &up->bcs, bias_list, eps_zero, kSq_mass, true, up->use_gpu); up->lhs_dual = gkyl_malloc(sizeof(double[up->numnodes_global])); @@ -477,6 +491,9 @@ void gkyl_fem_poisson_release(gkyl_fem_poisson *up) #endif } + if (up->bias_list_ho.num_bias_plane > 0) + gkyl_free(up->bias_list_ho.bp); + if (up->isdomperiodic) { gkyl_array_release(up->rhs_cellavg); gkyl_free(up->rhs_avg); diff --git a/moments/zero/gkyl_fem_poisson_priv.h b/moments/zero/gkyl_fem_poisson_priv.h index 176aca24f5..8284591317 100644 --- a/moments/zero/gkyl_fem_poisson_priv.h +++ b/moments/zero/gkyl_fem_poisson_priv.h @@ -1964,6 +1964,8 @@ struct gkyl_fem_poisson { // Objects to apply the LHS operator (M+K), used to apply 1-rho^2*Lap. bool has_lhs_apply; // Whether the LHS-apply objects are allocated. struct gkyl_poisson_bc bcs; // Saved BCs (used to build the mass solver). + struct gkyl_poisson_bias_plane_list bias_list_ho; // Host-side copy of the bias + // plane list (used to build the mass solver). struct gkyl_fem_poisson *mass; // Mass-matrix (M) solver for the LHS apply. double *lhs_dual; // Host global dual (M+K)*x vector. double *lhs_dual_cu; // Device global dual (M+K)*x vector. From 9c01440d058417dc45034a265e6daa5a4a929c0d Mon Sep 17 00:00:00 2001 From: Antoinehoff Date: Mon, 6 Jul 2026 10:13:30 -0700 Subject: [PATCH 08/12] activate flr effects in the tcv 2x2v adapt source regression test --- .../creg/rt_gk_tcv_iwl_adapt_source_2x2v_p1.c | 24 +++++++++++++------ 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/gyrokinetic/creg/rt_gk_tcv_iwl_adapt_source_2x2v_p1.c b/gyrokinetic/creg/rt_gk_tcv_iwl_adapt_source_2x2v_p1.c index 2e0ee0094a..809435d3e7 100644 --- a/gyrokinetic/creg/rt_gk_tcv_iwl_adapt_source_2x2v_p1.c +++ b/gyrokinetic/creg/rt_gk_tcv_iwl_adapt_source_2x2v_p1.c @@ -18,7 +18,7 @@ struct gk_app_ctx { // Geometry and magnetic field parameters double a_shift, Z_axis, R_axis, R0, a_mid, x_inner, r0, B0, kappa, delta, q0, Cy, Bref, x_LCFS; // Plasma parameters - double me, qe, mi, qi, n0, Te0, Ti0; + double me, qe, mi, qi, n0, Te0, Ti0, rho_i; // Collision parameters double nuFrac, nuElc, nuIon; // Source parameters @@ -353,6 +353,7 @@ struct gk_app_ctx create_ctx(void) double c_s = sqrt(Te0/mi); double omega_ci = fabs(qi*B0/mi); double rho_s = c_s/omega_ci; + double rho_i = vti/omega_ci; // Configuration domain parameters double Lx = Rmid_max-Rmid_min; // Domain size along x. @@ -446,6 +447,7 @@ struct gk_app_ctx create_ctx(void) .me = me, .qe = qe, .mi = mi, .qi = qi, .n0 = n0, .Te0 = Te0, .Ti0 = Ti0, + .rho_i = rho_i, .nuFrac = nuFrac, .nuElc = nuElc, .nuIon = nuIon, .num_sources = num_sources, @@ -625,6 +627,11 @@ main(int argc, char **argv) .collide_with = { "ion"}, }, + .flr = { + .type = GKYL_GK_FLR_PADE_CONST, + .Tperp = ctx.Te0, + }, + .source = { .source_id = GKYL_PROJ_SOURCE, .num_sources = ctx.num_sources, @@ -706,6 +713,11 @@ main(int argc, char **argv) .collide_with = { "elc"}, }, + .flr = { + .type = GKYL_GK_FLR_PADE_CONST, + .Tperp = ctx.Ti0, + }, + .source = { .source_id = GKYL_PROJ_SOURCE, .num_sources = ctx.num_sources, @@ -761,11 +773,6 @@ main(int argc, char **argv) .val = 0.0, // Biasing value. }, }; - - struct gkyl_poisson_bias_line_list bias_line_list = { - .num_bias_line = 2, - .bl = target_corner_bcs, - }; // Field. struct gkyl_gyrokinetic_field field = { @@ -775,8 +782,11 @@ main(int argc, char **argv) { .dir = 0, .edge = GKYL_LOWER_EDGE, .type = GKYL_BC_GK_FIELD_DIRICHLET, .value = {0.0} }, { .dir = 0, .edge = GKYL_UPPER_EDGE, .type = GKYL_BC_GK_FIELD_DIRICHLET, .value = {0.0} }, }, - .bias_line_list = &bias_line_list, .time_rate_diagnostics = true, + .flr = { + .type = GKYL_GK_FLR_PADE_CONST, + .gyroradius = ctx.rho_i, + }, }; // Geometry From 2ece5acdf3f23317beaf0ea2922fe3a58ef3dccb Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Mon, 6 Jul 2026 15:27:24 -0400 Subject: [PATCH 09/12] fix the flr species bc --- gyrokinetic/apps/gk_species.c | 22 ++++------------------ 1 file changed, 4 insertions(+), 18 deletions(-) diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index c59a4aedbd..612c89f68f 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/gyrokinetic/apps/gk_species.c @@ -1575,25 +1575,11 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st gks->flr_kSq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); gkyl_array_shiftc(gks->flr_kSq, -pow(sqrt(2.0),app->cdim), 0); // Sets kSq=-1. - // The gyroaverage uses varying-Dirichlet BCs (input field as its own - // boundary value) so no FLR correction is applied at the wall. - struct gkyl_poisson_bc flr_bc; + // Gyroaverage BCs: the input field is its own boundary value. + struct gkyl_poisson_bc flr_bc = app->field->poisson_bcs; for (int d=0; dcdim-1; d++) { - struct gkyl_gyrokinetic_bc *bc_lo = gk_fetch_bc_with_dir_edge(app->field->info.poisson_bcs, 2*app->cdim, d, GKYL_LOWER_EDGE); - if (bc_lo != 0) { - if (bc_lo->type == GKYL_BC_GK_FIELD_PERIODIC) - flr_bc.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); - else - flr_bc.lo_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_DIRICHLET_VARYING); - } - - struct gkyl_gyrokinetic_bc *bc_up = gk_fetch_bc_with_dir_edge(app->field->info.poisson_bcs, 2*app->cdim, d, GKYL_UPPER_EDGE); - if (bc_up != 0) { - if (bc_up->type == GKYL_BC_GK_FIELD_PERIODIC) - flr_bc.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_PERIODIC); - else - flr_bc.up_type[d] = gkyl_gyrokinetic_translate_poisson_bc_type(GKYL_BC_GK_FIELD_DIRICHLET_VARYING); - } + flr_bc.lo_type[d] = GKYL_POISSON_DIRICHLET_VARYING; + flr_bc.up_type[d] = GKYL_POISSON_DIRICHLET_VARYING; } // Deflated Poisson solve is performed on range assuming decomposition is *only* in z. gks->flr_op = gkyl_deflated_fem_poisson_new(app->grid, app->basis_on_dev, app->basis, From 07be7878821a28b178b8aee3978726ed66bc8e59 Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Mon, 6 Jul 2026 17:27:01 -0400 Subject: [PATCH 10/12] let the FLR BC be periodic if they need to be --- gyrokinetic/apps/gk_species.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index 612c89f68f..65d1172bd2 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/gyrokinetic/apps/gk_species.c @@ -1575,10 +1575,13 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st gks->flr_kSq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); gkyl_array_shiftc(gks->flr_kSq, -pow(sqrt(2.0),app->cdim), 0); // Sets kSq=-1. - // Gyroaverage BCs: the input field is its own boundary value. + // Gyroaverage BCs: the input field is its own boundary value at + // non-periodic boundaries. struct gkyl_poisson_bc flr_bc = app->field->poisson_bcs; for (int d=0; dcdim-1; d++) { + if (flr_bc.lo_type[d] != GKYL_POISSON_PERIODIC) flr_bc.lo_type[d] = GKYL_POISSON_DIRICHLET_VARYING; + if (flr_bc.up_type[d] != GKYL_POISSON_PERIODIC) flr_bc.up_type[d] = GKYL_POISSON_DIRICHLET_VARYING; } // Deflated Poisson solve is performed on range assuming decomposition is *only* in z. From b79e276e47365f5624d231bd58fd02694879d86f Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Fri, 24 Jul 2026 11:39:47 -0400 Subject: [PATCH 11/12] add a regression test for the FLR response --- gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c | 328 ++++++++++++++++++ 1 file changed, 328 insertions(+) create mode 100644 gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c diff --git a/gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c b/gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c new file mode 100644 index 0000000000..0cb0be5eae --- /dev/null +++ b/gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c @@ -0,0 +1,328 @@ +// Single-mode validation of the Pade FLR chain (see flr_effects.tex). +// +// Straight slab, uniform B, Dirichlet walls in x (the fully periodic +// perpendicular domain is a known unresolved issue in fem_poisson_perp), +// periodic in z, and a single ion density mode +// n_i = n0*(1 + A*sin(kx*x)), n_e = n0, kx = pi/Lx, +// which is an eigenfunction of every operator in the chain. The initial +// field solve is linear and diagonal in kx, so with b = (kx*rho_i)^2 the +// potential with FLR effects on and off satisfies +// phi_flr/phi_noflr = (1 + b)/(1 + b/2), +// from the density gyroaverage (Gamma_1 = 1/(1+b/2)) and the FLR +// inversion of phi (A = 1 + b). No time stepping is needed; this test +// applies ICs, solves the field, and checks amplitudes at frame 0. +// +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +struct flr_ctx { + int cdim, vdim; + double me, qe, mi, qi; + double n0, Te0, Ti0, B0; + double rho_i; // Ion gyroradius, sets the FLR scale. + double kx; // Perturbation wavenumber. + double pert_amp; // Relative amplitude of the density perturbation. + double Lx, Lz; + double vpar_max_elc, mu_max_elc, vpar_max_ion, mu_max_ion; + int cells[GKYL_MAX_DIM]; +}; + +static struct flr_ctx +create_ctx(void) +{ + double eV = GKYL_ELEMENTARY_CHARGE; + double me = GKYL_ELECTRON_MASS, mi = GKYL_PROTON_MASS; + double qe = -eV, qi = eV; + + double n0 = 1.0e19; + double Te0 = 100.0*eV, Ti0 = 100.0*eV; + double B0 = 1.0; + + double vte = sqrt(Te0/me), vti = sqrt(Ti0/mi); + double rho_i = sqrt(Ti0*mi)/(qi*B0); + + // One half-wavelength sine mode (vanishing at the x walls) with kx*rho_i = 0.5. + double kx_rho = 0.5; + double kx = kx_rho/rho_i; + double Lx = M_PI/kx; + double Lz = 1.0; + + struct flr_ctx ctx = { + .cdim = 2, .vdim = 2, + .me = me, .qe = qe, .mi = mi, .qi = qi, + .n0 = n0, .Te0 = Te0, .Ti0 = Ti0, .B0 = B0, + .rho_i = rho_i, + .kx = kx, + .pert_amp = 1.0e-2, + .Lx = Lx, .Lz = Lz, + .vpar_max_elc = 5.0*vte, + .mu_max_elc = me*pow(5.0*vte,2)/(2.0*B0), + .vpar_max_ion = 5.0*vti, + .mu_max_ion = mi*pow(5.0*vti,2)/(2.0*B0), + .cells = {32, 4, 12, 8}, + }; + return ctx; +} + +static void +density_elc(double t, const double *xn, double *fout, void *ctx) +{ + struct flr_ctx *app = ctx; + fout[0] = app->n0; +} + +static void +density_ion(double t, const double *xn, double *fout, void *ctx) +{ + struct flr_ctx *app = ctx; + double x = xn[0]; + fout[0] = app->n0*(1.0 + app->pert_amp*sin(app->kx*x)); +} + +static void +temp_elc(double t, const double *xn, double *fout, void *ctx) +{ + struct flr_ctx *app = ctx; + fout[0] = app->Te0; +} + +static void +temp_ion(double t, const double *xn, double *fout, void *ctx) +{ + struct flr_ctx *app = ctx; + fout[0] = app->Ti0; +} + +static void +zero_func(double t, const double *xn, double *fout, void *ctx) +{ + fout[0] = 0.0; +} + +static void +mapc2p(double t, const double *xc, double *xp, void *ctx) +{ + xp[0] = xc[0]; xp[1] = xc[1]; xp[2] = xc[2]; +} + +static void +bfield_func(double t, const double *xc, double *fout, void *ctx) +{ + struct flr_ctx *app = ctx; + fout[0] = 0.0; fout[1] = 0.0; fout[2] = app->B0; +} + +// Run the initial field solve with or without FLR effects and write frame 0. +static void +run_case(struct flr_ctx *ctx, struct gkyl_app_args *app_args, struct gkyl_comm *comm, + bool use_flr, const char *name) +{ + struct gkyl_gyrokinetic_species elc = { + .name = "elc", + .charge = ctx->qe, .mass = ctx->me, + .vdim = ctx->vdim, + .lower = { -ctx->vpar_max_elc, 0.0}, + .upper = { ctx->vpar_max_elc, ctx->mu_max_elc}, + .cells = { ctx->cells[2], ctx->cells[3] }, + .polarization_density = ctx->n0, + + .projection = { + .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM, + .ctx_density = ctx, + .ctx_upar = ctx, + .ctx_temp = ctx, + .density = density_elc, + .upar = zero_func, + .temp = temp_elc, + }, + + .bcs = { + { .dir = 0, .edge = GKYL_LOWER_EDGE, .type = GKYL_BC_GK_SPECIES_ABSORB, }, + { .dir = 0, .edge = GKYL_UPPER_EDGE, .type = GKYL_BC_GK_SPECIES_ABSORB, }, + }, + .num_diag_moments = 1, + .diag_moments = {GKYL_F_MOMENT_M0}, + }; + + struct gkyl_gyrokinetic_species ion = { + .name = "ion", + .charge = ctx->qi, .mass = ctx->mi, + .vdim = ctx->vdim, + .lower = { -ctx->vpar_max_ion, 0.0}, + .upper = { ctx->vpar_max_ion, ctx->mu_max_ion}, + .cells = { ctx->cells[2], ctx->cells[3] }, + .polarization_density = ctx->n0, + + .projection = { + .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM, + .ctx_density = ctx, + .ctx_upar = ctx, + .ctx_temp = ctx, + .density = density_ion, + .upar = zero_func, + .temp = temp_ion, + }, + + .bcs = { + { .dir = 0, .edge = GKYL_LOWER_EDGE, .type = GKYL_BC_GK_SPECIES_ABSORB, }, + { .dir = 0, .edge = GKYL_UPPER_EDGE, .type = GKYL_BC_GK_SPECIES_ABSORB, }, + }, + .num_diag_moments = 1, + .diag_moments = {GKYL_F_MOMENT_M0}, + }; + + if (use_flr) { + ion.flr.type = GKYL_GK_FLR_PADE_CONST; + ion.flr.Tperp = ctx->Ti0; + } + + struct gkyl_gyrokinetic_field field = { + .gkfield_id = GKYL_GK_FIELD_ES, + .polarization_bmag = ctx->B0, + .poisson_bcs = { + { .dir = 0, .edge = GKYL_LOWER_EDGE, .type = GKYL_BC_GK_FIELD_DIRICHLET, .value = {0.0} }, + { .dir = 0, .edge = GKYL_UPPER_EDGE, .type = GKYL_BC_GK_FIELD_DIRICHLET, .value = {0.0} }, + }, + }; + + if (use_flr) { + field.flr.type = GKYL_GK_FLR_PADE_CONST; + field.flr.gyroradius = ctx->rho_i; + } + + struct gkyl_gk gk = { + .cdim = ctx->cdim, + .lower = { 0.0, -ctx->Lz/2.0 }, + .upper = { ctx->Lx, ctx->Lz/2.0 }, + .cells = { ctx->cells[0], ctx->cells[1] }, + .poly_order = 1, + .basis_type = app_args->basis_type, + + .geometry = { + .geometry_id = GKYL_GEOMETRY_MAPC2P, + .world = {0.0}, + .mapc2p = mapc2p, + .c2p_ctx = ctx, + .bfield_func = bfield_func, + .bfield_ctx = ctx, + }, + + .num_periodic_dir = 1, + .periodic_dirs = {1}, + + .num_species = 2, + .species = { elc, ion }, + .field = field, + + .parallelism = { + .comm = comm, + .cuts = {app_args->cuts[0], app_args->cuts[1]}, + .use_gpu = app_args->use_gpu, + }, + }; + strcpy(gk.name, name); + + gkyl_gyrokinetic_app *app = gkyl_gyrokinetic_app_new(&gk); + gkyl_gyrokinetic_app_apply_ic(app, 0.0); + gkyl_gyrokinetic_app_write_field(app, 0.0, 0); + gkyl_gyrokinetic_app_release(app); +} + +// L2 norm of a DG field, sqrt(int f^2 dx). +static double +calc_l2(struct gkyl_rect_grid grid, struct gkyl_range range, struct gkyl_range range_ext, + struct gkyl_basis basis, struct gkyl_array *f) +{ + struct gkyl_array *l2 = gkyl_array_new(GKYL_DOUBLE, 1, range_ext.volume); + gkyl_dg_calc_l2_range(&basis, 0, l2, 0, f, range); + gkyl_array_scale_range(l2, grid.cellVolume, &range); + double l2red[1]; + gkyl_array_reduce_range(l2red, l2, GKYL_SUM, &range); + gkyl_array_release(l2); + return sqrt(l2red[0]); +} + +int +main(int argc, char **argv) +{ + struct gkyl_app_args app_args = parse_app_args(argc, argv); + struct flr_ctx ctx = create_ctx(); + + struct gkyl_comm *comm = gkyl_gyrokinetic_comms_new(app_args.use_mpi, app_args.use_gpu, stderr); + + run_case(&ctx, &app_args, comm, false, "rt_gk_flr_response_2x2v_p1_flroff"); + run_case(&ctx, &app_args, comm, true, "rt_gk_flr_response_2x2v_p1_flron"); + + // Read back both potentials. + struct gkyl_rect_grid grid; + gkyl_rect_grid_init(&grid, ctx.cdim, (double[]) { 0.0, -ctx.Lz/2.0 }, + (double[]) { ctx.Lx, ctx.Lz/2.0 }, (int[]) { ctx.cells[0], ctx.cells[1] }); + struct gkyl_basis basis; + gkyl_cart_modal_serendip(&basis, ctx.cdim, 1); + int nghost[GKYL_MAX_CDIM] = { 1, 1 }; + struct gkyl_range local, local_ext; + gkyl_create_grid_ranges(&grid, nghost, &local_ext, &local); + + struct gkyl_array *phi_off = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume); + struct gkyl_array *phi_on = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume); + struct gkyl_rect_grid grid_read; + gkyl_grid_sub_array_read(&grid_read, &local, phi_off, "rt_gk_flr_response_2x2v_p1_flroff-field_0.gkyl"); + gkyl_grid_sub_array_read(&grid_read, &local, phi_on, "rt_gk_flr_response_2x2v_p1_flron-field_0.gkyl"); + + // Measured amplitudes and expected factors. + double l2_off = calc_l2(grid, local, local_ext, basis, phi_off); + double l2_on = calc_l2(grid, local, local_ext, basis, phi_on); + + double b = pow(ctx.kx*ctx.rho_i, 2.0); + double ratio_expected = (1.0 + b)/(1.0 + b/2.0); + double ratio = l2_on/l2_off; + + // Analytic amplitude of the FLR-off solve: phi = qi*A*n0/(eps*kx^2)*cos(kx*x). + double eps_pol = ctx.n0*(ctx.mi + ctx.me)/pow(ctx.B0, 2.0); + double phi_amp = ctx.qi*ctx.pert_amp*ctx.n0/(eps_pol*pow(ctx.kx, 2.0)); + double l2_off_expected = phi_amp*sqrt(ctx.Lx*ctx.Lz/2.0); + + // Mode shape: phi_on must be the same mode scaled by the expected ratio. + struct gkyl_array *diff = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume); + gkyl_array_set(diff, 1.0, phi_on); + gkyl_array_accumulate(diff, -ratio_expected, phi_off); + double l2_shape = calc_l2(grid, local, local_ext, basis, diff); + + double err_ratio = fabs(ratio/ratio_expected - 1.0); + double err_abs = fabs(l2_off/l2_off_expected - 1.0); + double err_shape = l2_shape/l2_on; + + printf("\nFLR single-mode response test (kx*rho_i = %g, b = %g):\n", ctx.kx*ctx.rho_i, b); + printf(" |phi| FLR-off measured/expected: %e / %e (err = %.3e)\n", l2_off, l2_off_expected, err_abs); + printf(" |phi| FLR-on: %e\n", l2_on); + printf(" amplitude ratio measured/expected: %f / %f (err = %.3e)\n", ratio, ratio_expected, err_ratio); + printf(" mode-shape error: %.3e\n", err_shape); + + int num_fail = 0; + num_fail += err_ratio > 0.02; + num_fail += err_abs > 0.05; + num_fail += err_shape > 0.02; + printf("%s\n", num_fail == 0? "PASSED" : "FAILED"); + + gkyl_array_release(phi_off); + gkyl_array_release(phi_on); + gkyl_array_release(diff); + gkyl_gyrokinetic_comms_release(comm); + + return num_fail; +} From 7c9e0b4e274062d8734c86dd73c22c8bc6e2070a Mon Sep 17 00:00:00 2001 From: antoinehoff Date: Fri, 7 Aug 2026 15:37:10 -0400 Subject: [PATCH 12/12] update the flr reg test to verify rho_i<<1 and rho_i>>1 limits, this passes. I have to come up with a less trivial case now. --- gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c | 51 ++++++++++++++----- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c b/gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c index 0cb0be5eae..e3c07ba48e 100644 --- a/gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c +++ b/gyrokinetic/creg/rt_gk_flr_response_2x2v_p1.c @@ -12,6 +12,11 @@ // inversion of phi (A = 1 + b). No time stepping is needed; this test // applies ICs, solves the field, and checks amplitudes at frame 0. // +// We sweep kx*rho_i (and thus b) +// from small to large to check two limits of the Pade chain: +// - b -> 0: rho_i -> 0 limit, ratio -> 1 (FLR effects vanish). +// - b >> 1: ratio -> 2, the large-argument saturation of the Pade form. +// #include #include #include @@ -43,7 +48,7 @@ struct flr_ctx { }; static struct flr_ctx -create_ctx(void) +create_ctx(double kx_rho) { double eV = GKYL_ELEMENTARY_CHARGE; double me = GKYL_ELECTRON_MASS, mi = GKYL_PROTON_MASS; @@ -56,8 +61,6 @@ create_ctx(void) double vte = sqrt(Te0/me), vti = sqrt(Ti0/mi); double rho_i = sqrt(Ti0*mi)/(qi*B0); - // One half-wavelength sine mode (vanishing at the x walls) with kx*rho_i = 0.5. - double kx_rho = 0.5; double kx = kx_rho/rho_i; double Lx = M_PI/kx; double Lz = 1.0; @@ -257,16 +260,17 @@ calc_l2(struct gkyl_rect_grid grid, struct gkyl_range range, struct gkyl_range r return sqrt(l2red[0]); } -int -main(int argc, char **argv) +static int +run_and_check_case(int case_idx, double kx_rho, struct gkyl_app_args *app_args, struct gkyl_comm *comm) { - struct gkyl_app_args app_args = parse_app_args(argc, argv); - struct flr_ctx ctx = create_ctx(); + struct flr_ctx ctx = create_ctx(kx_rho); - struct gkyl_comm *comm = gkyl_gyrokinetic_comms_new(app_args.use_mpi, app_args.use_gpu, stderr); + char name_off[64], name_on[64]; + snprintf(name_off, sizeof(name_off), "rt_gk_flr_response_2x2v_p1_case%d_flroff", case_idx); + snprintf(name_on, sizeof(name_on), "rt_gk_flr_response_2x2v_p1_case%d_flron", case_idx); - run_case(&ctx, &app_args, comm, false, "rt_gk_flr_response_2x2v_p1_flroff"); - run_case(&ctx, &app_args, comm, true, "rt_gk_flr_response_2x2v_p1_flron"); + run_case(&ctx, app_args, comm, false, name_off); + run_case(&ctx, app_args, comm, true, name_on); // Read back both potentials. struct gkyl_rect_grid grid; @@ -281,8 +285,11 @@ main(int argc, char **argv) struct gkyl_array *phi_off = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume); struct gkyl_array *phi_on = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume); struct gkyl_rect_grid grid_read; - gkyl_grid_sub_array_read(&grid_read, &local, phi_off, "rt_gk_flr_response_2x2v_p1_flroff-field_0.gkyl"); - gkyl_grid_sub_array_read(&grid_read, &local, phi_on, "rt_gk_flr_response_2x2v_p1_flron-field_0.gkyl"); + char fname_off[96], fname_on[96]; + snprintf(fname_off, sizeof(fname_off), "%s-field_0.gkyl", name_off); + snprintf(fname_on, sizeof(fname_on), "%s-field_0.gkyl", name_on); + gkyl_grid_sub_array_read(&grid_read, &local, phi_off, fname_off); + gkyl_grid_sub_array_read(&grid_read, &local, phi_on, fname_on); // Measured amplitudes and expected factors. double l2_off = calc_l2(grid, local, local_ext, basis, phi_off); @@ -322,6 +329,26 @@ main(int argc, char **argv) gkyl_array_release(phi_off); gkyl_array_release(phi_on); gkyl_array_release(diff); + + return num_fail; +} + +int +main(int argc, char **argv) +{ + struct gkyl_app_args app_args = parse_app_args(argc, argv); + struct gkyl_comm *comm = gkyl_gyrokinetic_comms_new(app_args.use_mpi, app_args.use_gpu, stderr); + + // Sweep kx*rho_i from the rho_i->0 limit (ratio -> 1) to the large-b Pade saturation (ratio -> 2). + double kx_rho_list[] = { 0.05, 0.1, 0.25, 0.5, 1.0, 2.0, 4.0 }; + int num_cases = sizeof(kx_rho_list)/sizeof(kx_rho_list[0]); + + int num_fail = 0; + for (int i = 0; i < num_cases; i++) + num_fail += run_and_check_case(i, kx_rho_list[i], &app_args, comm); + + printf("\n%d/%d cases failed.\n", num_fail, num_cases); + gkyl_gyrokinetic_comms_release(comm); return num_fail;