It appears that parallel simulations, with or without shared memory, are not working when using NonUniformRectCart.

I'm seeing errors like

Fatal error in PMPI_Group_incl: Invalid rank, error stack:
PMPI_Group_incl(184).............: MPI_Group_incl(group=0x88000000, n=4, ranks=0x402cb4c8, new_group=0x402cb5e8) failed
MPIR_Group_check_valid_ranks(256): Duplicate ranks in rank array at index 1, has value 0 which is also the value at index 0

Hi @pcagas, I'm currently working on implementing a plasma sheath model in a five moment multi-fluid code. In this paper and your dissertation, a vacuum boundary condition is imposed in the ghost cells for the two-fluid model. Could you please point me towards which Riemann Solver is used for these boundary conditions?

Vlasov objects are growing in complexity. I propose that we separate them, especially in the app system.

For example, when you look at the MaxwellField app there are lots of if-statements where it checks whether the plasmaB exists, if it does then it assumes VM, otherwise it assumes VP. Maybe we should just have separate Maxwell and Poisson apps. Just an idea.

Of course there are other models that add further complexity (SR Vlasov, PKPM) in both the field and the species code. Some of this could be simplified with separation.

Reason for this proposal: to reduce risk of breaking one model when modifying another.

So far our approach to making BGK collisions conservative in GK and VM has been to use a Lagrange fix (lagFix) so that the moments of the distribution function that comes out of MaxwellianOnBasis (MOB) matches the moments inputed into MOB. There are two problems with this approach:

1. It can change the distribution function in an undesirable way, e.g. increases it near the velocity boundaries.
2. The 3x2v lagFix kernel takes forever to generate (> 18 hours), indicating that it'll be impractically big.

There may be other issues with the lagFix implementation for GK/GkBGK also because 1x2v p=2 GkBGK appeared to diverge when using lagFix.

Not sure how to proceed here.

a. We could implement a series of rescalings in GkBGK in the spirit of what is done for fixing the initial conditions in GkProjection.
b. If we are willing to accept the distortion of the lagFix, we could consider implementing them in new kernels that inver the linear system with Eigen instead of Maxima. This may lead to a smaller 3x2v kernel.

Create shared libraries that are then linked to the final executable. At present, all object files are being linked to the final binary, making it huge and bloated.

Also allow specification of a path where user-define shared objects files are located, so these can be loaded by gkeyll.

Also remove p=2 primMom kernels for 3x3v.

This will make the build easier/faster on some systems.

Hello,

We recently encountered an issue with running Gkeyll on the newest version of Arch Linux. If Gkeyll is initialised using OpenMPI and more than 2 workers, then the simulation will fail. Or, if 1 to 2 workers are used on a non-uniform grid then the grid file is not written to a local file.

For those looking for a fix, we downgraded packages (all packages on our local machines) to their state at 30/07/2023 which solves the issue.

Presumably this is an error relating to "adios 1.13.1.post5" interacting with the new glibc library (2.38-1 and later). A file containing the full output is below.

error_output.txt

The primary point here is that we need a way to pass function pointers into __global__ CUDA kernels. these function pointers will point to particular existing kernels (instrumented with __host__ __device__ CUDA preprocessors). this way we can have generic __global__ kernels that do loops over cells, and then in each cell we can call a particular kernel function.

The DistFuncMomentCalc updater is the furthest along on GPU implementation, and should be used as a template in many regards. Credit to @manauref. So I will use this as an example for what I mean. I'm going to go through the currently implemented DistFuncMomentCalc workflow step by step and then make some comments.

Here is the CUDA-relevant block in Updater/DistFuncMomentCalc.lua where things start:

   if GKYL_HAVE_CUDA and self.calcOnDevice then
      local d_PhaseGrid  = grid:copyHostToDevice()
      local d_PhaseRange = Range.copyHostToDevice(phaseRange)

      local deviceNumber     = cudaRunTime.GetDevice()
      local deviceProps, err = cudaRunTime.GetDeviceProperties(deviceNumber)

      local phaseRangeDecomp = LinearDecomp.LinearDecompRange {
         range = phaseRange:selectFirst(pDim), numSplit = grid:numSharedProcs() }
      local numCellsLocal = phaseRange:selectFirst(pDim):volume()

      local numThreads = GKYL_DEFAULT_NUM_THREADS
      local numBlocks  = math.floor(numCellsLocal/numThreads)+1

      self._momCalcFun(d_PhaseGrid, d_PhaseRange, deviceProps, numBlocks, numThreads, distf:deviceDataPointer(), mom1:deviceDataPointer())

      cudaRunTime.Free(d_PhaseRange)
   else
      .... regular code ....

The C function that is pointed to by self._momCalcFun is cuda_MomentCalc1x1vSer_M0_P1. It lives in Updater/momentCalcData/DistFuncMomentCalcSer1x1vDevice.cu. This is a C function that wraps a CUDA kernel call:

// C function that wraps call to moment calculation global CUDA kernel
void cuda_MomentCalc1x1vSer_M0_P1(RectCart_t *grid, Range_t *range, GkDeviceProp *prop, int numBlocks, int numThreads, const double *fIn, double *out) {

  int warpSize = prop->warpSize;

  d_calcMom1x1vSer_M0_P1<<<numBlocks, numThreads, 2*(numThreads/warpSize)*sizeof(double)>>>(grid, range, fIn, out);
}

The global CUDA kernel d_calcMom1x1vSer_M0_P1 also lives in Updater/momentCalcData/DistFuncMomentCalcSer1x1vDevice.cu. It threads over cells, as we have discussed before. Each thread is responsible for a cell, and it calls MomentCalc1x1vSer_M0_P1. This is one of the typical maxima-generated C kernels that we know and love, except that it has been instrumented with __host__ __device preprocessors to allow it to be called from within this __global__ CUDA kernel.

This workflow works. However, as is it would require separate reduction kernels and wrapper kernels for every moment, dimensionality, polynomial order, etc. So this is where we need to do a bit more work to generalize these functions so that they can be reused. (I'm sure that Mana has thought about these things already, but just wanted to write up a full example.)

I think the global kernel which calls the local kernel and does the reduction should be generic so that it can be used for any moment, dimensionality, polynomial order, etc. Such a generic kernel would then need some way to know what local kernel function to call (in this case MomentCalc1x1vSer_M0_P1). Perhaps templating is the way to go for this, but I'm not sure if CUDA kernels can be templated. In any case, I'm thinking about something like this:

• The lua function pointer self._momCalcFun could point to a wrapper to some generic moment reduction kernel, something like

// C function that wraps call to moment calculation global CUDA reduction kernel
void cuda_calcMoment_wrapper(RectCart_t *grid, Range_t *range, GkDeviceProp *prop, int cdim, int pdim, function_pointer, int numBlocks, int numThreads, const double *fIn, double *out) {

  int warpSize = prop->warpSize;

  cuda_calcMoment<<<numBlocks, numThreads, 2*(numThreads/warpSize)*sizeof(double)>>>(grid, range, cdim, pdim, function_pointer, fIn, out);
}

where we now are directly passing cdim and pdim, and also a pointer to a moment calculator kernel function function_pointer, which in this specific case would point to MomentCalc1x1vSer_M0_P1.

• The global kernel __global__ cuda_CalcMoment will be more or less the same as what Mana has in d_calcMom1x1vSer_M0_P1, except for the hard-coded dimensional parameters will be replaced with parameters related to cdim and pdim, and the local kernel function (currently hard-coded to be MomentCalc1x1vSer_M0_P1) will be taken from the function_pointer argument. This will allow this global kernel to be reused for every moment, dimensionality, polynomial order, etc.

This structure will be similar for other updaters, where we will want to have a generic looping kernel, into which we pass dimensionality, polynomial order, etc, and also a function pointer to the particular local kernel function (our maxima-generated kernels) that should be called. We just need to figure out the architecture for passing function pointers into kernels.

At present some C++ headers have Gk prefix, some Gkyl and some none at all. We need to unify the prefix and naming scheme. It is very confusing otherwise.

function MomentSpecies:updateInDirection(dir, tCurr, dt, fIn, fOut, tryInv)
   local status, dtSuggested = true, GKYL_MAX_DOUBLE
   local tryInv_next = false
   if self.evolve then
      self:applyBc(tCurr, fIn, dir)
      assert(self:checkInv(fIn))

Line 495 of MomentSpecies.lua raises an error in... some cases, but because there's no descriptive text on the error and no documentation on what fIn or checkInv, it's impossible to tell what exactly causes the error. I'm guessing there's some issue with either a non-invertible matrix somewhere, or with the definition of equationInv in the species definition.

Heres my code that's causing the error. (Please note -- this is a bit of a mishmash of different examples at the moment, and because I'm still learning, some of the initial conditions may be nonphysical.)

https://pastebin.pl/view/42736129

When run, the console prints

*** LOAD ERROR ***
 ...e/andrea/gkylsoft/gkyl/bin/App/Species/MomentSpecies.lua:495: assertion failed!

with no further information.

Dear Gkeyll developers,

I am a researcher in plasma wakefield acceleration, one of the developers of the particle-in-cell (PIC) code WarpX https://ecp-warpx.github.io, and we are looking for an open-source plasma simulation code for several problems of interest to us, and have a few questions to see if Gkeyll could be the way to go. While PIC works like a charm to simulate plasma acceleration below a picosecond, we are currently interested in longer-term plasma evolution (up to a microsecond), where other models (in particular two-fluid) could provide a good description.

The image below shows what we are interested in: a high-energy particle beam perturbs an initially uniform neutral plasma, and drives an electron plasma wave (the "wake") in which gigantic accelerating fields can be used to accelerate another particle beam located in the first bubble behind the driver (the rightmost electron bucket in blue). We are currently interested in how this plasma wave relaxes long after the driver's passage. The important features for this process would be:

• The problem is axisymmetric, so an RZ description would be ideal.
• Collisions are likely to play a role long after the driver's passage, in particular electron-ion collisions.
• The background plasma is pre-ionized, with an ionization fraction of 1-10%, so collisional ionization is probably important.

Additional comments:

• The figure is in normalised units. In SI units, the beam would be ~10 microns wide, ~50 microns long, and the plasma density would be on the order 1.e15-1.e17 cm-3.
• Just modelling Hydrogen would be extremely helpful.
• The actual setup has some transverse boundaries, but simulating an infinitely-wide plasma would be a great first step.

Do you think Gkeyll could help for our problem?

Thanks in advance!

A high-energy, high-current electron beam (red) propagates (from left to right) at the speed of light through a neutral plasma and perturbs the electron density (blue), driving a non-linear wake. Ions are immobile on this short time scale, and start moving later on.

This discussion could also be interesting to @GregoryBoyle @SeverinDiederichs @MathisMewes.

I followed the directions, but ran into an error:

Before I try to debug it, are there any quick fixes already known?

git clone [email protected]:ammarhakim/gkyl.git
cd gkyl/
./machines/mkdeps.adroit.sh
./machines/configure.adroit.sh
./waf build install

Waf: Entering directory `/home/henryfs/gkyl/build'
[  1/782] Compiling Lib/lfs.c
[  2/782] Compiling Lib/base64.cpp
[  3/782] Compiling Lib/whereami.c
[  4/782] Compiling Lib/gkyl_ipow.cpp
Waf: Leaving directory `/home/henryfs/gkyl/build'
Build failed
Traceback (most recent call last):
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Task.py", line 180, in process
    ret=self.run()
  File "<string>", line 27, in f
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Task.py", line 173, in exec_command
    return self.generator.bld.exec_command(cmd,**kw)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Context.py", line 181, in exec_command
    ret,out,err=Utils.run_process(cmd,kw,cargs)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Utils.py", line 598, in run_process
    return run_prefork_process(cmd,kwargs,cargs)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Utils.py", line 545, in run_prefork_process
    raise OSError(trace)
WafError: Execution failure: ['icpc', '-O3', '-std=c++17', '-fPIC', '-ILib', '-I../Lib', '-ICuda', '-I../Cuda', '-I../../gkylsoft/luajit/include/luajit-2.1', '../Lib/gkyl_ipow.cpp', '-c', '-o/home/henryfs/gkyl/build/Lib/gkyl_ipow.cpp.1.o']
Traceback (most recent call last):
  File "<string>", line 30, in run
  File "/usr/lib64/python2.7/subprocess.py", line 711, in __init__
    errread, errwrite)
  File "/usr/lib64/python2.7/subprocess.py", line 1327, in _execute_child
    raise child_exception
OSError: [Errno 2] No such file or directory


Traceback (most recent call last):
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Task.py", line 180, in process
    ret=self.run()
  File "<string>", line 27, in f
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Task.py", line 173, in exec_command
    return self.generator.bld.exec_command(cmd,**kw)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Context.py", line 181, in exec_command
    ret,out,err=Utils.run_process(cmd,kw,cargs)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Utils.py", line 598, in run_process
    return run_prefork_process(cmd,kwargs,cargs)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Utils.py", line 545, in run_prefork_process
    raise OSError(trace)
WafError: Execution failure: ['icpc', '-O3', '-std=c++17', '-fPIC', '-ILib', '-I../Lib', '-ICuda', '-I../Cuda', '-I../../gkylsoft/luajit/include/luajit-2.1', '../Lib/base64.cpp', '-c', '-o/home/henryfs/gkyl/build/Lib/base64.cpp.1.o']
Traceback (most recent call last):
  File "<string>", line 30, in run
  File "/usr/lib64/python2.7/subprocess.py", line 711, in __init__
    errread, errwrite)
  File "/usr/lib64/python2.7/subprocess.py", line 1327, in _execute_child
    raise child_exception
OSError: [Errno 2] No such file or directory


Traceback (most recent call last):
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Task.py", line 180, in process
    ret=self.run()
  File "<string>", line 27, in f
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Task.py", line 173, in exec_command
    return self.generator.bld.exec_command(cmd,**kw)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Context.py", line 181, in exec_command
    ret,out,err=Utils.run_process(cmd,kw,cargs)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Utils.py", line 598, in run_process
    return run_prefork_process(cmd,kwargs,cargs)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Utils.py", line 545, in run_prefork_process
    raise OSError(trace)
WafError: Execution failure: ['icc', '-O3', '-Wall', '-fPIC', '-ILib', '-I../Lib', '-ICuda', '-I../Cuda', '-I../../gkylsoft/luajit/include/luajit-2.1', '../Lib/lfs.c', '-c', '-o/home/henryfs/gkyl/build/Lib/lfs.c.1.o']
Traceback (most recent call last):
  File "<string>", line 30, in run
  File "/usr/lib64/python2.7/subprocess.py", line 711, in __init__
    errread, errwrite)
  File "/usr/lib64/python2.7/subprocess.py", line 1327, in _execute_child
    raise child_exception
OSError: [Errno 2] No such file or directory


Traceback (most recent call last):
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Task.py", line 180, in process
    ret=self.run()
  File "<string>", line 27, in f
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Task.py", line 173, in exec_command
    return self.generator.bld.exec_command(cmd,**kw)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Context.py", line 181, in exec_command
    ret,out,err=Utils.run_process(cmd,kw,cargs)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Utils.py", line 598, in run_process
    return run_prefork_process(cmd,kwargs,cargs)
  File "/home/henryfs/gkyl/.waf-2.0.19-1f3c580272b15a03d2566843c5fe872a/waflib/Utils.py", line 545, in run_prefork_process
    raise OSError(trace)
WafError: Execution failure: ['icc', '-O3', '-Wall', '-fPIC', '-ILib', '-I../Lib', '-ICuda', '-I../Cuda', '-I../../gkylsoft/luajit/include/luajit-2.1', '../Lib/whereami.c', '-c', '-o/home/henryfs/gkyl/build/Lib/whereami.c.1.o']
Traceback (most recent call last):
  File "<string>", line 30, in run
  File "/usr/lib64/python2.7/subprocess.py", line 711, in __init__
    errread, errwrite)
  File "/usr/lib64/python2.7/subprocess.py", line 1327, in _execute_child
    raise child_exception
OSError: [Errno 2] No such file or directory

Also, have you considered CMake instead of custom scripts for every machine? If you need help getting started, I can point you at a book about using Modern CMake (that I wrote).

(Okay, I promise this issue isn't me missing a note or having a corrupted compiler this time. I know because I was able to work around it.)

The following error concerns the pre-g0 branch.

When attempting to fetch ADIOS 1.13.1 during mkdeps, curl instead fetches an html error file. gunzip of course has trouble interpreting this, and the build of ADIOS halts.

gzip: adios-1.13.1.tar.gz: not in gzip format
tar: adios-1.13.1.tar: Cannot open: No such file or directory
tar: Error is not recoverable: exiting now
./build-adios.sh: line 14: cd: adios-1.13.1: No such file or directory
./build-adios.sh: line 15: ./configure: No such file or directory
make: *** No rule to make target 'install'.  Stop.

This is because the current url used to fetch ADIOS 1 is has broken SSL certs. (See line 11 of build-adios.sh)

curl -L https://users.nccs.gov/~pnorbert/adios-1.13.1.tar.gz > adios-1.13.1.tar.gz

It can be downloaded manually through the browser by changing "http" to "https". However, that's not enough to satisfy curl though, which still raises an error if you change the url.

curl: (60) SSL certificate problem: unable to get local issuer certificate
More details here: https://curl.se/docs/sslcerts.html

curl failed to verify the legitimacy of the server and therefore could not
establish a secure connection to it. To learn more about this situation and
how to fix it, please visit the web page mentioned above.

This could probably be fixed by linking to a different mirror of ADIOS, or changing the security requirements on curl. It would likely be a relatively simple fix, just not one I'm experienced enough in net security to do myself.

Workaround:

• Delete or comment out lines 9 and 11 of build-adios.sh (the lines starting with rm and curl), to prevent the script from replacing files you add manually.
• Download ADIOS from https://users.nccs.gov/~pnorbert/adios-1.13.1.tar.gz and place it in your install-deps folder.
• Run ./install-deps/mkdeps.sh or ./machines/mkdeps.[OS].sh as appropriate.

it throws the error

GkProjection.lua:58: attempt to index field 'species' (a nil value)

I think because fullInit, or some other method, that should be called for the global population is only being called for the local population. That, or the comm/rank options in the :read method are interfering with this.

We are nearly ready for renaming g0-merge to main and using that for all simulations. However, we are not yet there yet. Here is a list of milestone we should meet before we do this rename. Please comment/discuss.

• Remove all unused C++ kernels and code from the system
• Cleanup usage of ZeroArray so we do not need to check for GPUs every time we make a call to the G0 code. This is very confusing and I do not think it is really needed.
• Cleanup and simplify CartField. Many features are not needed and so remove them.
• Switch to Git version of LuaJIT so we are on the latest LJ code and not some older, forked version.
• Remove dependencies that are no longer needed for G2. For example, is Eigen needed? I not, lets remove. (The dispersion Tool use Eigen but I think this is easy to replace with OpenBLAS functions).
• Re-enable automatic regression tests to ensure we have not broken anything. Even if this is on local machines we need to do this before we rename.

A secondary thing we can do but perhaps can wait is:

• Wrap Moments code in using Lua-C API (no FFI but straight C). We will need this for fluid sims very soon.

At the moment comparefiles works like

gkyl comparefiles -a fileA.bp -b fileB.bp

However, for the purpose of comparing data serially generated by our regression tests (i.e. accepted results) with data generated by a parallel run it'd be good to be able to do something like

gkyl comparefiles -da directoryA/ -db directoryB/

and it'll compare files just like our regression system does, reporting relative differences. This is a useful step on the way to checking results with MPI.

A more specific example, suppose I generate data for the vm-lbo/rt-lboRelax-1x1v-p1 test with
gkyl runregression run create -r vm-lbo/rt-lboRelax-1x1v-p1.lua
and make some changes that primarily affect MPI, then i'd like to check if a parallel run produces the same results with

mpirun -np N gkyl vm-lbo/rt-lboRelax-1x1v-p1.lua
gkyl comparefiles -da $HOME/gkylsoft/gkyl-results/vm-lbo/rt-lboRelax-1x1v-p1/ -db vm-lbo/

mpirun -n 1 gkyl -g -e "tFinal=2.e-6;numFrames=8;wallTime=172800." input-file.lua

produces

*** LOAD ERROR ***
[string "GKYL_USE_GPU = truetFinal=2.e-6;numFrames=8;w..."]:1: unexpected symbol near '='

Updater/Bc is old and should be removed, especially since it doesn't run on GPUs.

We need to replace its use in MaxwellField.

In the g0-merge branch, when running on GPUs, ASheathPotential is using Mpi.Bcast to broadcast the sheath potential. I think this either implicitly relies on CUDA-aware MPI or is just broadcasting the wrong potential. It should be revised.

See fe841e6

Basically, we currently have two ways of importing sources from a file, via fromFile or profile. The former does not scale by the gyro-center coordinate Jacobian, while the latter does. This ought to be addressed so source importing is more consistent.

I'm thinking the source written out should be always be written without the Jacobian. Then both options (fromFile and profile) would be made to always scale by the Jacobian.

Oftentimes, a Vlasov simulation might fail due to too low of collisionality. In these cases the last output is usually garbage (including the restart frame for all Dynvector quantities) and not what you want to restart from, so you rename a frame that looks reasonable to have the suffix restart, and then restart the simulation at higher collisionality.

This worked before when all Dynvector quantities (and cflFrac) were written out as separate files, but now that they are a single file it does not work.

I am not immediately sure what the solution should be. Perhaps an option to generate completely new dynvector files and leave it to the user to merge them together in post-processing. I am currently testing if I just comment out the read of the dynVector restart file if I can at least restart and take time-steps, but a more elegant solution is likely desired.

I am still trying to figure out if the last saved time-step size is also an issue, as the *dt.bp file appears to be empty. ReadRestart is already doing a dummy forwardEuler call. Perhaps a solution is to have an option to just get the size of the stable time-step from this dummy forwardEuler call if no *dt.bp file is detected

Some GkLBO tests may be performing suboptimally. For example, the gk-lbo/rt-lboRelax-1x2v-p1 test is taking nearly 8x longer than vm-lbo/rt-lboRelax-1x2v-p1. It would be worthwhile understanding the origin of this.

1x2v simulations run in parallel on NERSC's Cori with one MPI process per cell are throwing the error

*** LOAD ERROR ***
...a/cori/gkeyll/gkylsoft/gkyl/bin/DataStruct/CartField.lua:138: cannot convert 'nil' to 'double'

when requesting the following boundary flux diagnostics

diagnosticBoundaryFluxMoments = {"GkM0", "GkM1", "GkM2", "GkUpar", "GkTemp", "GkEnergy"},
diagnosticIntegratedBoundaryFluxMoments = {"intM0", "intM1", "intKE", "intHE"},

Weirdly I can't reproduce it on my laptop.

I've tracked it down to

self._confBoundaryFieldPtr[c] = inFldPtr[c]

in the :evalOnConfBoundary method of Updater/Bc.lua.

I don't yet know the solution, but I have noticed one potential memory problem. The :evalOnConfBoundary method creates a CartField inside an isFirst if-statement. However, for diagnostics GkBeta, GkEnergy, GkUparCross, GkVtSqCross I believe :evalOnConfBoundary gets called again before isFirst gets flipped (which happens in the :advance method). I interpret this to mean that a CartField gets created where a CartField already existed.

There is an inconsistency in the way CartField handles non-numeric types. In particular:

• Sync() is not done correctly when the data stored in the field is non-numeric. At present, non-numeric types are sent as byte-arrays but the sizeof(elem) is not used, which means data is lost.

• It is assumed in the accumulate etc methods that the field stores doubles.

• It is assumed when doing a copy to GPU that the field only stores doubles. This is potentially serious as it means we can't use floats or any other number types via the CartField on the GPU. Should change this. (It may not really matter at present as all kernels assume doubles anyway).

We need to refactor the field object. We need to separate the data store and the operations on the data. This will simplify the CartField code, which at present is very complex (though < 1000 lines of code!)

When running with useShared = true, if the decompCuts does not equal the number of nodes, one still obtains the error that:

ERROR: CartProdDecomp: Number of sub-domains and processors must match

but the simulation does not time out and will still keep running. Likely due to issue that shared communicator thinks it's still valid and tries to share memory across compute nodes.

Error present as of commit c15e2df on March 31st. Likely still present as shared memory has had issues lately. Need a lot more testing of shared memory.

Description of the bug

In the g0-merge, ZeroArray.Array.__gc() is never called.

Impact

This is causing memory leak, but only at the end of a simulation when all the bulk data arrays go out of scope. In other word, it perhaps doesn't affect simulation results.

Possible causes

The cdata being wrapped, struct gkyl_array is created by an external function call gkyl_array_new, not by ffi itself through new. Therefore, no implicit ffi.gc call is automatically done during the creation of an instance. See this post.

Possible solutions

• Solution 1, Create the instance through ffi.new, then call some init function to initialize it.
  • This new function gkyl_array_init can use codes of the main body of present gkylzero.gkyl_array_new. The latter can then call gkyl_array_init to avoid duplicate codes.
  • How to handle the gpu copy then?
  • This could be generalized to other gkylzero classes, if needed. That is, we may always prefer to have an independent _init method that initializes a given object.
• Solution 2, following Ammar's workaround, by wrapping over the following struct gkyl_array_container rather than directly struct gkyl_array. Then, in the instance creation, allocate a new gkyl_aray_container object using ffi.new, and embed the actual gkyl_array created by gkyl_array_new.
  • This requires wider code changes, including all of ZeroArray member methods, and any code that uses CartField._zero itself (i.e., a struct gkyl_array), needs to be written as CartField._zero.arr.
  • See 5650513 and 07a1c91

struct gkyl_array_container {
   struct gkyl_array *arr;
};

local array_fn = {
   copy = function (self, src)
      return ffiC.gkyl_array_copy(self.arr, src.arr)
   end,
...
}

local array_mt = {
   __new = function (self, atype, ncomp, size, on_gpu)
      local container = new(self)
      container.arr = ffiC.gkyl_array_new(atype, ncomp, size)
      return container
   end,
   __gc = function(self)
      self:release()
   end,
   __index = array_fn,
}

@nmandell @ammarhakim @JunoRavin

I disabled test6 in Unit/test_CartFieldCuda.lua because it was throwing an error in Portal (see below). Perhaps this is not an issue and I'm running the test the wrong way. If that's the case, can the author enlighten me?

Also, this unit test prints out a bunch of indices. Is this left over from some dev or debug work? should it be removed?

To allow multiblock geometries and also allow (initially) static mesh refinement (i.e. put patches of finer grid on coarse grids) we need GridCollection and FieldCollection objects. With this, we need to modify updaters to loop over each grid in GridCollection and update the equation on that grid. Some subtle issue will remain with moments etc.

For this, a refactoring of the CartField object is needed. It is horridly complex at present and needs to be simplified. Moving the MPI stuff out of it will help.

(disclaimer, my understanding of fluid sources is incomplete)
(note, given the input from a weekend slacker I made some edits in the first comment below)

Currently sources, in the sense of App/Sources and the sources table in PlasmaOnCartGrid, are only used by fluids. They apply CollisionlessEmSource and TenMomentRelaxSource (and Jimmy/Liang have a TenMomentGradSource in a the gradient-closure branch). Comments/observations:

1. ``sources`` are currently separated from the rest of the fluid update (``updateInDirection``), and applied (in ``updateSources`` with some fractional time step) before and after (or is it just before/after when using super-time-stepping?) ``updateInDirection`` (in the ``fvDimSplit`` time stepper).
2. In an input file, these fluid sources appear as tables outside of species tables. Even if the source is just a closure term that applies to a single species. This is just because we want to operator-split these terms and treat them as a ``source``.
3. I find it counter intuitive to call all these things sources. I see you are calling them that way because if you write the equations in conservative form, they appear as terms on the RHS (actually this isn't even true for the heat flux in the energy equation). This is unlike what in tokamak fluid and gyrokinetic turbulence codes people call sources (in my experience), which are typically terms that correspond to plasma/energy/momentum sources into the system. I propose we "rename" some, especially as ``sources`` is starting to be used for other even less source-like terms (like closures) just because of the functional location of the ``sources`` table and the operator splitting happening in the time stepper. This is not too dissimilar from something Jimmy or Liang considered, to rename the folder App/Sources to something like App/FluidXXX.

An aside: I'm not a big fan of packing multiple Apps in the same file (e.g. GkField/GkGeometry, MaxwellField/ExternalMaxwellField).

I would like to propose we refactor sources (or several Apps for that matter) using something along these lines:

1. Create something like an App/Forces folder (pls. propose other names too). In that folder we could place CollisionlessEmSource (rename Lorentz Force?) and ExternalMaxwellField.
2. TenMomentRelaxSource, TenMomentGradSource and Diffusion could go into a App/Dissipation folder (pls. propose other names too).
3. Give these Apps, and the collisions Apps as well, another input e.g. applyOp = "pre", "post", "both", "none"/nil. These options indicate:
   • "pre": the operator only gets applied at the beginning of a time step, before the other RK stages, explicit terms, or the rest of the fluid update.
   • "post": like "pre" but applied at the end instead.
   • "both": operator is applied at the beginning and at the end, e.g. with dt/2 time step each, or we could also request (or code default) pre and post dtFracs if necessary.
   • "none"/nil: no operator splitting takes place, and these terms are applied (presumably explicitly) inside the species:advance().
     
     Not all operators will have all of these options, e.g. CollisionlessEmSource may only support applyOp="both".
4. In input files:
   4.1 single-species operators (ones that get computed independently for each species) like Diffusion, Collisions (TenMomentGradSource?) could live inside the species tables.
   4.2 Operators that get computed for all species at once (CollisionlessEmSource) could either
   • live in each species table (with some duplication of code in input files),
   • go in the field table (after all, this App is just dictating how the fields get applied, and one could in the future envision a situation in VM in which the fields are applied implicitly, thus requiring something like CollisionlessEmSource),
   • or we could leave it outside the species tables.
     
     Then we could sort out what gets applied when in PlasmaOnCartGrid, in a similar way to how collisions get sorted out during kinetic species initialization. That is, we would look through species and field tables looking for "sub-Apps" with the input applyOp. If applyOp="pre"/"post"/"both", we put them in preOps and/or postOps tables accordingly. App with applyOp="none"/nil get placed in a corresponding collisions/sources/etc tables inside each species.
5. In fvDimSplit, instead of looping over sources before and after updateInDirection, we would instead loop over preOps and postOps. And inside species we would only advance/apply the collisions/sources/etc that had applyOp="none"/nil.
6. In forwardEuler or the RK steppers, we could also loop over preOps and postOps for cases in which we are applying, for example, diffusion or collisions implicitly/super-time-stepped.

Some motivations:

• Modify App/Sources so that we can put some kinetic sources in there.
• Make changes to PlasmaOnCartGrid to bring fluids and kinetics a little closer.
• Implement support for implicit/supper-time-stepped treatment of diffusion/collisions in RK time steppers.

I have cleaned up the rununit command of the runregression Tool. To use this do:

gkyl runregression rununit

From the Regression directory.

Please add this to the nightly regression system. To do this do:

gkyl runregression rununit > outUnit
cat outUnit | grep FAILED!
cat outUnit | grep FAILED! | wc -l

Which will run the unit tests, list the failed unit tests and count the number of failed tests. Please add this information to the nightly regression emails.

We have been looping over the species tables with loops like

for nm, s in pairs(species) do
...
end

pairs is the correct iterator for a table of name-value pairs. However, it is not deterministic: there is no guarantee that consecutive runs will loop in the same order. More importantly, there's no guarantee that different MPI ranks will loop in the same order. This is a problem for I/O (and possibly other operations), as one rank may be trying to output the ion distribution function at the same time as another rank may be outputting the electron distribution function.

I can conceive of 3 options for fixing this:

A) Create a list of sorted keys, which we can store in the species table:

local species_keys = {}
for k in pairs(species) do table.insert(species_keys, k) end
table.sort(species_keys)
species["keys"] = lume.clone(species_keys)

then we can loop like

for i = 1, #species["keys"] do
   local nm, s = species["keys"][i], species[species["keys"][i]]
   ...
end

B) Store the key table as the metadata table of species. First create species_keys as in option A above. Then

setmetatable(species,species_keys)

then loop like

for _, nm in ipairs(getmetatable(species)) do
   local s = species[nm]
   ...
end

C) Create a centrally located iterator. For example we could create the speciesIter iterator and place it in gky/Lib, so that we would simply loop with

local speciesIter = require "Lib.speciesIter"

for nm, s in speciesIter(species) do
   ...
end

Option A has been implemented in branch newSpeciesLoops, commit 6ec2195. I'm waiting to merge to get input from others. The more I think about it, I'm not a fan of option C and if we are not going to use the metatable property of species for anything else, I'd vote for option B.

This issue apparently arose when we switched from global LF fluxes, to local, quadrature point-based LF fluxes.

Dear All,

Thanks for providing the code to the public.
I received error as (after first running the code):

(base) [usr@localhost vm-damp]$ gkyl vm-damp.lua
Tue Apr 27 2021 10:49:46.000000000
Gkyl built with
Gkyl built on Apr 24 2021 18:42:00
Initializing Vlasov-Maxwell simulation ...
** WARNING: timeStepper not specified, assuming rk3
ERROR: Attribute element changeset has invalid value attribute
ERROR: Attribute element changeset has invalid value attribute
[localhost:04570] *** Process received signal ***
[localhost:04570] Signal: Segmentation fault (11)
[localhost:04570] Signal code: Address not mapped (1)
[localhost:04570] Failing at address: 0xfffffff7
[localhost:04570] [ 0] /lib64/libpthread.so.0(+0xf370)[0x7f0a7c6c0370]
[localhost:04570] [ 1] /lib64/libc.so.6(cfree+0x1c)[0x7f0a7b3a638c]
[localhost:04570] [ 2] gkyl(adios_common_define_attribute_byvalue+0x3ce)[0x42f248]
[localhost:04570] [ 3] gkyl(adios_define_attribute_byvalue+0x50)[0x4290be]

I would be thankful to help me to correct run the first example.

All The Best

Regression tests rt-vm-init.lua and rt-vm-init-p2.lua have an issue with the Lagrange Fix, seemingly due to the driftSpeed being a vector function in higher than 1V. When running the following error is produced:
ProjectOnBasis: this updater was created for a scalar/vector field, not a vector/scalar field.

This issue is to document the problem, but we may be deprecating the LagrangeFix updater owing to it not always being useful because of how it changes the tails of the distribution to fix the energy. If iterative schemes for conservative BGK work, issue may be closed as we will be able to officially deprecate LagrangeFix

I finally figured out what the issue with the code I linked in #174 was, and it's caused by another, entirely separate, mostly-silent error.

The parameter polyOrder is unused in the Moments App, but the app doesn't ignore it or raise an error if it's defined. Instead, it tries to accommodate both the polynomial order you gave it and its own solving methods. This leads to NaNs infesting your higher order moments, and the same cryptic "assertion failed" error I made an issue about a bit ago. Even assuming that that error had a more descriptive message ("Invalid pressure value: NaN"), it would still lead you to think that it was an issue with your initial conditions, and not with a hanging parameter.

A simple check to compare the App type with whether or not this parameter was defined would help a lot -- whether it made sure that polyOrder's value didn't affect the simulation, or if it just raised an error when it was defined. Either would work.

I've been migrating code from a Vlasov-Maxwell setup to the 10-moment app, and I must have missed that one line in the process. The only way I was able to figure out that that was the cause was by copying whole chunks of a working app over piece by piece.

Currently Gkeyll has an issue with how M1 Macs deal with denormalized floating point numbers. There is some key functionality that is removed due to this.

Currently, 4 unit tests fail. test_TenMomenSrc.lua, test_VlasovRect.lua, test_EvalOnNodes.lua, test_CartGrid.lua.

To build, in gkyl.cxx in the main directory, one must comment out - _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON); and - fesetenv(FE_DFL_DISABLE_SSE_DENORMS_ENV);

In Updater/wscript, delete the following lines

FemPerpPoissonImpl.cpp
FemParPoissonImpl.cpp
FemMatrices.cpp
FemGyroaverageImpl.cpp
EvalOnNodesImpl.cpp

Running the vlasov quickstart parallel example from here:

https://gkeyll.readthedocs.io/en/latest/quickstart/vlasovExample1.html#running-the-simulation

using the lua script here:

https://gkeyll.readthedocs.io/en/latest/quickstart/inputFiles/vm-tsw-2x2v.html

on a workstation using mpich 3.3.2 fails with the following error:

$ mpirun -np 10 gkyl vm-tsw-2x2v.lua                                                                                      
Thu Oct 06 2022 15:15:25.000000000
Gkyl built with 4c3e5680196b
Gkyl built on Sep 22 2022 12:37:06
Initializing Vlasov-Maxwell simulation ...
*** LOAD ERROR ***
 /lore/cwsmith/spaceWeather/gkylsoft/bin/Comm/Mpi.lua:389: bad argument #3 to 'MPI_Win_shared_query' (cannot convert 'unsigned int [1]' to 'int64_t *')
 <... snip ...>

On an Intel Mac we also tried letting the gkyl install scripts install OpenMPI and hit the same errors.

On the workstation, making the following changes to disable shared memory and increase the number of subdomains appeared to avoid the failures:

101,102c101,102
<    decompCuts = {1,1},                      -- Cuts in each configuration direction
<    useShared = true,                        -- If using shared memory
---
>    decompCuts = {5,2},                      -- Cuts in each configuration direction
>    useShared = false,                       -- If using shared memory

My setup is able to run VlasovMaxwell and Gyrokinetic apps just fine.

However, trying to run any of the given example scripts for Euler/Ten-moment sims gives the following error:

> gkyl rt-10m-gem.lua
*** LOAD ERROR ***
 [string "-- Gkyl -------------------------------------..."]:2: attempt to call field 'Moments' (a nil value)

I set up this small test script.

local Apps = require("App.PlasmaOnCartGrid")

for key,value in pairs(Apps) do
    print(key, value)
end

Running it through gkyl gives the following:

> gkyl test.lua
Gyrokinetic     function: 0x7fe173393770
VlasovMaxwell   function: 0x7fe1733937e0

There's no entry for a Moment app! Somehow, my version of Gkeyll only has VlasovMaxwell and Gyrokinetic apps available. Redownloading and rebuilding Gkeyll didn't help.

Found this out while testing some p=0 code. Presumably EvalOnNodes in this case should just evaluate at cell center.

I'm running sims with this field solver in the input file:

   -- Field solver.
   field = Plasma.Field {
      epsilon0 = epsilon_0, mu0 = mu_0,
      evolve   = true, -- Evolve field?
      hasMagneticField = false,
      bcLowerPhi = {{T="D", V=V_left}},
      bcUpperPhi = {{T="D", V=0.0}},
   },

This should output a 4 component output with phi and 3 components of 0s (since no B field being included). However, all the outputs are the electric potential.

Redo some of the Maxima code to reduce duplication in GK, Vlasov, Maxwell etc equations. Likely this needs a set of common functions to generate code for hyperbolic equations and also some refactoring of recovery code.

Dear Gkeyll team,

Pardon my somewhat vague question. We are in the starting phase of simulating kinetic effects relating to instabilities in the ionosphere, and consider using Gkeyll's Vlasov solver for this. More specifically, we'd like to pass a beam through a plasma experiencing one of the instabilities well known to appear in the auroral region. There is a range of different instabilities that could be of interest, so we could start with one that is easier to simulate using Gkeyll. Some examples could be cyclotron instabilities, current convective instabilities, inhomogeneous energy density-driven instability (IEDDI), Kelvin-Helmhotz instability (KHI) and gradient drift instabilities (GDI) (although I expect these latter to be difficult to simulate due to the larger scale). My question is whether any of you have any experience to suggest which instability would be easier to simulate using Gkeyll?

How many cells do we need?

• In 1d, a five-cell stencil is needed to update each real cell.
• This is because two ghost cells are required on each side to update the real cell with 2nd-order accuracy.
• More exactly, the two ghost cells on one side are used to compute waves at the ghost-ghost cell face and real-ghost cell face. The former is only used to the latter, and the latter is needed to update the adjacent real cell.

One way to do cell-based update

• If we do a pure cell-based update (update one real cell per thread) with numCellsPerBlock=numThreadsPerBlock,
  • on the block, the number of Riemann problems to be solved is numRp=numCellsPerBlock+3.
  • One possibility is to request the threads with threadIdx.x==0 and threadIdx.x==blockDim.x to do the three additional Rp to get the additional waves needed.
  • The computed waves are stored in shared memory for computing the smoothness parameter for flux limiting.

Another way

• We may also require each block to load a chunk of the data with two ghost cells on each side.
• Then each thread solves a Riemann problem to compute all waves needed, including those used for computing smoothness only.
  • In this case, numExtendedEdgesPerBlock=numThreadsPerBlock.
• This way is like decompose the local domain on the cpu onto even smaller domains on the blocks.
• This way, an alternative invIndex mapping is needed (and might be very complicated).

Several tools no longer work. In particular, the tools needing linear algebra libraries are not working. These should be fixed. Also, seems man/woman keyword searches are somewhat wonky (for example gkyl man maxima gives the wrong web-page). Fix this also.

We are ready to replace master with g0-merge, at least on the Vlasov and gyrokinetic side.

Dear Gkeyll team, first of all, thanks for making this software.

I am trying to build gkyl on a new cluster, but get the following eight errors (and many warnings, of which I show three) when running waf (in verbose mode):

[1203/1308] Compiling Eq/vlasovData/VlasovStreamMax1x2vP2.cpp
12:18:32 runner ['/cluster/software/CUDAcore/11.1.1/bin/nvcc', '-x', 'cu', '-c', '-dc', '-O3', '-std=c++14', '-lineinfo', '--maxrregcount=255', '-arch=sm_70', '--compiler-options="-fPIC"', '-IEq', '-I../Eq', '-ICuda', '-I../Cuda', '-IGrid', '-I../Grid', '-ILib', '-I../Lib', '-IBasis', '-I../Basis', '-IEq', '-I../Eq', '-IDataStruct', '-I../DataStruct', '-IEq/vlasovData', '-I../Eq/vlasovData', '-IEq/maxwellData', '-I../Eq/maxwellData', '../Eq/vlasovData/VlasovStreamMax1x2vP2.cpp', '-o', '/cluster/projects/nn9299k/software/gkyl-2.0_source/build/Eq/vlasovData/VlasovStreamMax1x2vP2.cpp.2.o']
../Eq/vlasovData/VlasovTmplModDecl.h(564): error: identifier "VlasovVolStream3x3vMaxP2" is undefined

../Eq/vlasovData/VlasovTmplModDecl.h(574): error: identifier "VlasovSurfStream3x3vMax_X_P2" is undefined

../Eq/vlasovData/VlasovTmplModDecl.h(577): error: identifier "VlasovSurfStream3x3vMax_Y_P2" is undefined

../Eq/vlasovData/VlasovTmplModDecl.h(580): error: identifier "VlasovSurfStream3x3vMax_Z_P2" is undefined

../Eq/vlasovData/VlasovTmplModDecl.h(586): error: identifier "VlasovVol3x3vMaxP2" is undefined

../Eq/vlasovData/VlasovTmplModDecl.h(597): error: identifier "VlasovSurfElcMag3x3vMax_VX_P2" is undefined

../Eq/vlasovData/VlasovTmplModDecl.h(600): error: identifier "VlasovSurfElcMag3x3vMax_VY_P2" is undefined

../Eq/vlasovData/VlasovTmplModDecl.h(603): error: identifier "VlasovSurfElcMag3x3vMax_VZ_P2" is undefined

8 errors detected in the compilation of "../Eq/vlasovData/VlasovTmplModDeviceFuncs.cpp".

../Eq/vlasovData/VlasovSurfElcMagMax3x3vP1.cpp(14): warning: variable "dv1" was declared but never referenced

../Eq/vlasovData/VlasovSurfElcMagMax3x3vP1.cpp(14): warning: variable "wv1" was declared but never referenced

../Eq/vlasovData/VlasovSurfElcMagMax3x3vP1.cpp(17): warning: variable "B0" was declared but never referenced

gkyl worked well on my laptop not long ago so I realize I'm probably doing something wrong. Do you have any idea what might cause these errors, or how I can go about debugging it?

For reference, I'm using commit 9267eeff17dee3be03f80816d1ad7c2004ed7f12. I should also mention that I have changed the name of the gkylsoft folder in the machine files to better match the naming convention our group uses on the cluster.

Thank you.

I'm attempting to install Gkeyll for Ubuntu via WSL. My system is a HP Pavilion Aero Laptop 13-be1xxx, running Windows 11.

I'm able to set up dependencies and run the configuration just fine. When I attempt to compile Gkeyll itself (./waf build install), I get the following error:

[16/24] Creating build/gkylgittip.h
[18/24] Linking build/gkyl
/usr/bin/ld: /home/andrea/gkylsoft/gkylzero/lib/libgkylzero.so: undefined reference to `LAPACKE_dgesv'
collect2: error: ld returned 1 exit status

Waf: Leaving directory `/home/andrea/Sources/gkyl/build'
Build failed
 -> task in 'gkyl' failed with exit status 1 (run with -v to display more information)

I've tried deleting and reinstalling dependencies, as well as building from source instead of from machine files. I always run into this error.

Am I missing some package that I need to apt-get? Is there a config option I'm missing?

I think there's an erroneous memory access in line 71 of Grid/MappedCart.lua.

It is assigning a value to d1[3] and d2[3], but these vectors have size 2, not 3. I suspect that the correct code that should replace lines 69-71 is:

   d1[1], d2[1], d3[1] = diff.derivt(self._mapc2p, 1)(xc)
   d1[2], d2[2], d3[2] = diff.derivt(self._mapc2p, 2)(xc)

Otherwise the current code takes the derivative with respect to a third variable that doesn't exist. If this is correct, I wonder why test_2() in Unit/test_MappedCart.lua works. Maybe I'm wrong and this issue should be resolved/closed.

Allow restarts to happen on arbitrary number of MPI ranks.

The dataset inside a DynVector sometimes change shape after a restart. Noah tried to prepare postgkyl to handle this (postgkyl commit 89dbc7d5936bd31a4a26e0857bd5261f63cd996f), but it doesn't seem to always work.

This leads to errors like

pgkyl -f gk24-wham1x2v_esEnergy.bp pl
...
ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)

and inspecting the file with bpls we find

bpls gk24-wham1x2v_esEnergy.bp
...
  double   TimeMesh713   {2292}
  double   Data713       {2292}
  double   TimeMesh714   {2292}
  double   Data714       {2292, 1}

Also, I've also seen the insertion of empty datasets like:

bpls gk24-wham1x2v_elc_intM0.bp
...
  double   TimeMesh1019  {0}
  double   Data1019      {0, 1}
  double   TimeMesh1020  {1}
  double   Data1020      {1, 1}
  double   TimeMesh1021  {1}
  double   Data1021      {1, 1}
  double   TimeMesh1022  {0}
  double   Data1022      {0, 1}

but I'm not sure if this is related.