Refactor phi wall boundary conditions into species table - #1100
Closed
Maxwell-Rosen wants to merge 68 commits into
Closed
Refactor phi wall boundary conditions into species table#1100Maxwell-Rosen wants to merge 68 commits into
Maxwell-Rosen wants to merge 68 commits into
Conversation
…e maximum magnetic field is determined. Bmag max is stored as a gkyl_array. Right now, we only do this for bmag, but we need to store phi as a 1d maximum array. I haven't decided on the final design for how 2x OAP simulations should be accomodated. Perhpas we need some general gkyl_dg_array_reduce methods that take a 2D array and turn it into a 1D array instead of a 0D number. This is a kind of reduction method, but it's not a total reduction. I tested the regression test included for a 2x2v boltzmann mirror and the output of the magnetic field looks correct. The current implementation evaluates bmag at cell corners, but we ideally should do the corners in Z, but the quadrature nodes in psi.
- Introduced a new header file `gkyl_array_dg_find_peaks.h` that defines a structure and functions for finding peaks (local maxima, minima, and boundary values) in DG fields. - Implemented an internal structure in `gkyl_array_dg_find_peaks_priv.h` to manage peak finding operations, including storage for peak values and coordinates. - Removed unused initialization and writing of `bmag_max` arrays in `gyrokinetic.c` to streamline the geometry setup process. - Deleted the `gkyl_gk_geometry_bmag_max_init` and `gkyl_gk_geometry_bmag_max_release` functions from `gk_geometry.c` as they are no longer needed, simplifying the geometry management.
…es another DG array at the peaks of the initilized array. Add appropriate unit tests which pass to ctest_array_dg_find_peaks. Update gk_species_fdot_multiplier to use the project_on_peaks function with phi. Now, everything passed to loss_cone_mask_gyrokinetic is a gkyl_array. The loss_cone_mask is updated accordingly
…ount of compution we need for evaluating phi at its peak in the app. Unit tests pass. Regression tests look fine as well. They're all valgrind clean. I think the right way to do the paralellism is to do the peak finding on a global bmag, just like how it is done for the position_map, then when we evaluate phi, all processes evaluate it at this peak, however only one will return a true value. This process will broadcast the array to the rest of the processes
…on at the peak locations
… method, which is just like find_peaks, but it computes the global maximum or minimum.
…nd regression tests are brought over from another branch. Unit tests for the array mask, loss cone mask, and the regression tests for the kinetic electron POA mirror are valgrind free
…grind clean. Regression test is added and produces reasonable results.
…formance - Introduced a helper function `mkarr` to streamline array allocation for GPU and CPU. - Removed the `gkyl_loss_cone_mask_gyrokinetic_Dbmag_quad_wall` function and integrated its logic into the main processing flow. - Updated the GPU kernel `gkyl_loss_cone_mask_gyrokinetic_Dbmag_quad_cu_ker` to compute `Dbmag_quad` directly from `bmag_peak` instead of `bmag_max`. - Enhanced tandem mirror support by adding handling for `bmag_peak` and `phi_m` in the GPU kernels. - Simplified the logic for determining trapped particles in the `gkyl_loss_cone_mask_gyrokinetic_ker` and `gkyl_loss_cone_mask_gyrokinetic_quad_ker` functions. - Improved readability and maintainability by restructuring conditional checks and variable assignments.
…plier. Add possibility of kinetic electrons and tandem mirrors. The damping regression test is failing, both here and on main. They are for different issues. Main fails because the loss_cone updater has an issue. Here, it fails because it's using scale_by_cell with a multi-component array. I'm not sure the right way to fix this
…ove the aspects about the cellwise evaluation and quadrature points because that breaks the array_scale_by_cell method which is used
… refactor this in the future, but it is just proof of concept for now to make sure it works correctly.
… arrays need to be passed to objects like the loss_cone_mask, where it expects these to be GPU arrays. It's just easier to have this module fit the archatecture of the rest of the code, rather than doing something different and copying between device and host. It wouldn't interface well. Claude generated most of the cuda code, with strong guidence from Maxwell
…h compute sanitizer with the array_find_peaks which was causing crashes in the loss cone mask. These issues are fixed. There was some funny business regarding the basis being on the host vs device. Refactor the allocations in the GPU kernels to not be inside the kernels. Instead, it's allocated at init time. The GPU code pulses, which is odd, but it runs. The 2x2v POA regression test runs and is compute sanitizer clean. The other POA tests do not error either on GPU and are compute sanitizer clean.
…itting code to main and broke a unit test with the geometry enum changes
…hi_smooth_global array for improved performance and consistency across computations.
…ion of doing the kinetic electron and tandem mirror. The code is built again to make sure nothing is affected.
…=2 relevant code for the peak finders
…egression test (and in my production simulations)
…and switch to user-defined damping type in simulation context.
…d enhance comments for clarity.
- Removed unused includes and commented-out code for clarity. - Simplified the initialization of quadrature values by consolidating the logic into a single function. - Introduced a new function to compute escape barriers based on the potential and magnetic field. - Replaced the previous approach of handling quadrature nodes with a more efficient nodal representation. - Updated the main advance function to utilize the new escape barrier calculations and streamline the trapped particle detection logic. - Cleaned up memory management by ensuring proper release of allocated arrays.
…remove unnecessary parameters
…nce communication logic for global data assembly
… parallel execution with MPI
…and streamline node handling logic
- Removed unused inline functions for node coordinate matching and field value calculations in loss_cone_mask_gyrokinetic.c. - Simplified the escape_barriers function by directly integrating its logic into the main kernel. - Updated the CUDA kernel for loss cone mask gyrokinetic to streamline the computation of escape barriers and Hamiltonian checks. - Enhanced the handling of configuration and phase ranges in the CUDA implementation for better clarity and performance. - Adjusted the kernel launch parameters to accommodate the new structure of the code and ensure proper execution.
…nd enhance mask handling
…ul, but let's not clutter the codebase
…-run allgather Root cause The reset path itself was fine â�� the bug was the NCCL allgather of bmag being re-issued mid-run. The loss-cone init in gk_species_fdot_multiplier.c did: gkyl_comm_array_allgather(app->comm, &app->local, &app->global, app->gk_geom->geo_corn.bmag, fdmul->bmag_global); bmag is static geometry, but gk_species_fdot_multiplier_reset does a full release â�� init, so this collective fired again on every reset. On Perlmutter's new stack (cudatoolkit/13.0 + nccl/2.29.2-cu13 + AWS Libfabric/GDR), issuing that collective from a reset corrupts the run â�� matching every observation: - 1 GPU fine (no real collective) / 4 GPU all NaN / 8 GPU edge ranks NaN / both phases right after a reset / didn't exist before the CUDA13/NCCL upgrade. - Without reset it's fine because the only mid-run collective that breaks is the reset's re-gather â�� the per-step phi allgather and the app-creation bmag allgather both work. On the host/device question: I checked, and geometry is already on device at app-creation fdot init too â�� app_new runs app_new_geom (hostâ��device swap at gyrokinetic.c:442-446) before app_new_solver creates species. So host-vs-device wasn't the difference; the new NCCL just doesn't tolerate that mid-run gather. The fix Gather the static global bmag once (in the reliable app-construction context) and reuse it across resets, so no allgather is ever issued from a reset: - struct gk_species gains fdot_bmag_global (cached, gathered lazily on first loss-cone init). - Loss-cone (re)inits now gkyl_array_acquire that cache instead of re-gathering. - Released at species teardown; refcounting keeps reset release/init balanced. After this, no fdot reset issues any collective (FIXED_FACTOR never did; loss-cone now reuses the cache) â�� which is exactly why it should clear the "both phases, right after reset" failure. 3 files, +21/â��5; CPU build is clean (0 errors). I can't run the NCCL path here, so please rebuild with configure.perlmutter.gpu.sh and re-test the 4-GPU case (the cleanest repro: it went from all-NaN to expected to be fine). One thing to watch: if any of your real input decks first enable the loss-cone multiplier via a mid-run reset (rather than having it on at app creation like failure.c), the one-time gather would then happen mid-run. If you want it bullet-proof for that case too, say the word and I'll make the cache gather eagerly during species construction regardless of the initial multiplier type
…void mid-run allgather" This reverts commit 9c6d202.
- Implemented a new gyrokinetic simulation code for a mirror configuration that supports two species: ions and electrons. - Introduced a pseudo orbit-averaged (POA) integrator with distinct phases for orbit averaging and full dynamics. - Defined plasma parameters, geometry, and source terms for the simulation context. - Added functions for magnetic field calculations, density evaluations, and velocity mappings for both species. - Included diagnostics for integrated quantities and phase-space outputs. - Established a framework for running multiple phases with adjustable parameters and time-stepping controls.
…n on multiple nodes / GPUs on the latest NCCL and CUDA. It has nothing to do with NCCL or CUDA!!!! it's becuase time_dilation_scale_const can be initilized to nan when using gkyl_malloc, which happens when the poa_phases are initilized, however the regression tests do not pass this parameter into the poa_phases struct, so the regression tests were fine, but my production simulations were not. using calloc initilizes everything to zero, which leads to an elegant handling of this parameter. Here is claude's final response. It was never NCCL â�� it was uninitialized memory, and MPI just exposed it. In rt_gk_failure.c, poa_phases is gkyl_malloc'd (plain malloc, not zeroed) and the phase loop never sets time_dilation_scale_const. So run_phase feeds an uninitialized value into the loss-cone multiplier: - Serial / 1-GPU â�� poa_phases lands on fresh OS pages (zeroed) â�� field is 0.0 â�� the old if (x <= 0.0) guard catches it â�� 1.0 â�� runs fine (your = 1 runs). - MPI â�� MPI_Init dirties the heap first, so the field reuses garbage holding a NaN bit pattern â�� -nan. -nan <= 0.0 is false, so it slips through, multiplies into the multiplier, and dt becomes -0/NaN. Each rank gets different garbage â�� looked like a "some regions / comms" bug for the whole investigation. Fixes applied (both compile clean) 1. Guard in gk_species_fdot_multiplier.c (what you asked for): falls back to 1.0 for any non-finite or non-positive time_dilation_scale_const, with a one-time warning. - Crucial subtlety: the build uses -ffast-math, so isfinite()/x!=x get optimized away (clang even warned -Wnan-infinity-disabled). The guard therefore detects non-finite via the bit pattern ((u >> 52) & 0x7ff == 0x7ff), which fast-math can't elide. 2. Root cause in rt_gk_failure.c: gkyl_malloc â�� gkyl_calloc for poa_phases, so unset fields (time_dilation_scale_const, f_threshold, the _ion/_elc variants) are zero instead of garbage. You need to rebuild â�� your three runs above were the old binary. With the rebuild, the 4-GPU run should print the warning (or nothing, with the calloc fix) and proceed with scale = 1.0 instead of NaN. One thing worth deciding: with time_dilation_scale_const correctly defaulting to 1.0, the OAP gets no time dilation. If the production sim intends a specific dilation factor during OAP, set it explicitly in the phase loop (poa_phases[2*i].time_dilation_scale_const = ...) â�� otherwise the calloc/guard just makes "unset" mean "1.0" deterministically. Want me to wire an explicit OAP value through, or leave it at 1.0?
…tiplier
- Introduced new enum `gkyl_loss_cone_boundary_type` to classify boundary behaviors (open, sheath, closed).
- Updated `gk_species_fdot_multiplier` to handle lower and upper wall potentials for sheath boundaries.
- Modified `gkyl_loss_cone_mask_gyrokinetic` to incorporate wall potentials in escape barrier calculations.
- Enhanced tests to cover various boundary condition scenarios, ensuring correct behavior for both open and closed boundaries.
- Adjusted CUDA kernel to accommodate new parameters for wall potentials.
- Sheath barriers now use the distinct material-wall potential, including grounded,
biased, asymmetric, and evolving walls.
- Corrected path scanning to include both cell traces and domain endpoints.
- Added roundoff-safe cutoff equality handling.
- Wired species BCs and global wall potentials through gyrokinetic/apps/
…ement for species
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This is from sol 5.6 high. The parameter for the shetath potential function should be a property of the species, not the field, says Mana. This is because the sheath BC is a species BC, not a field BC.
This is behind my PR #997
Sol's notes