Test this function. The difficulty here is that one has to make a hierarchy with at least one refined level and advance the levels so the time interpolation can be checked.
The test should check the ghost nodes are ok at different sub-steps of the refined level to check it's doing the time interpolation correctly.

AdvanceLevel

This issue is about implementing the SolverPPC::avdanceLevel.
This method advances the hybrid state from time n to time n+1 via the predictor-predictor-corrector scheme.

The method is decomposed in three main steps (red) each calling sub-steps:

• predictor1

  • faraday (blue) updates the magnetic field
  • ohm (green) updates the electric field
  • average fields at half step (yellow)

• predictor2

  • faraday updates B
  • ohm updates E
  • average E and B at half step

• corrector

  • faraday updates B
  • Ohm updates E

After each predictor, there's a call to the IonUpdater to update the ions particles (only after predictor 2) and moments (after both predictors).

hybrid initialize method does not compute ion moments yet.
As a result the root level moments are not initialized

• call core::computeIonMoments in initialize()
• check moments are computed OK by comparing to the user inputs

We should not compute the resistive or hyper-resistive term if their coef is zero.
And we should not test the coef at run time while evaluating the ohm's law.

There should code that plugs these terms only if needed.

One idea maybe is to put say each term is a Functor, inheriting from an abstract functor "OhmsLawTerm". The object Ohm would have a collection of these abstract terms and computing ohm's law could be something like:

(simplfied version...) :

for (ix=psi: ix <=pei; ++ix)
{
    for (auto& term : terms)
    {
        Ex(ix) += term(ix);
    }
}

therefore upon building the Ohm object, we would only add the relevant terms in the equation, discarding resistive and hyper-resistive terms if their coef is zero.

• Push my commits to GitHub
• Finish my changes
• Open a pull request

Objective

The objective of this issue is to implement the Electrons.
Electrons is an data object, like Ions, it is part of the HybridState.
It is used in the hybrid PPC solver to

• get the electron density
• get the bulk velocity
• get electron pressure

Electron moments could a priori be obtained from different assumptions. The most frequent are:

• quasi-neutrality, which says that:

  • Ne = Ni
  • Ve = Vi - J/(Ne*e)

• isothermal pressure closure

  • which says that the electron temperature Te is a spatio-temporal constant
  • Pe = Ne*Te

The pressure could be obtained from another closure equation, such as adiabatic, polytropic, etc. More rarely, the electron density and bulk velocity could also be obtained from other assumptions (like evolving the electron fluid equation, or evolving particles etc.).

The was moments are calculated should therefore, a priori, not be known at the solver level.

Most probably, the different possible ways of getting electron moments will need informations from other HybridState quantities.

Problem ?

In the PPC solver, starting from all quantities known at time n, electronMomentModel is called after Faraday and Ampère, just before Ohm. Hence, we have Ne = Ni at n, Vi at n and J at n+1. Then, Ve = Vi-J/N is calculated with a mix of quantities given at either n, or n+1. It should not be a big issue as it is needed for the calculation of E in Ohm, but it is not very nice. ..

Find and implement a good test for HybridHybridMessengerStrategy::synchronize().
A priori it concerns:

• make sure the algos and schedules work as expected to sync data
• make sure schedules still work (are remade) after regriding

parameters to change for tests:

• change number of levels
• change number of refine boxes per level
• change the position of the refine boxes for each level, including corner cases
• loop over all interpolation orders
• loop over all dimensions (for later)
• serial VS MPI run

initialization tests:

These tests are performed after the initialization of a full Hybrid Hierarchy.

Electromagnetic initialization

• L_0 E-M fields are EXPECT_DOUBLE_EQ with the user supplied functions on all
  physical nodes. user supplied functions are different for each field and non constant
  user supplied function are periodic in case of periodic boundary condition.

• in case of periodic boundary conditions : L0 patch ghost on the boundary of
  the domain are equal to their domain counterpart

• patch ghost nodes are equal to their interior counterpart on neighbor L_0 patches

• L_i (i!=0) E-M fields are equal to the user supplied function in case of a linear
  function (that is different for each field) on physical nodes

• L_i (i!=0) E-M fields patch ghosts nodes are EXPECT_DOUBLE_EQ to their counterpart
  in neighbor L_i patches (and not spatially constant)

• L_i (i!=0) E-M fields level ghost nodes are filled with values equal to the interpolation
  of the underlying coarser level (equal to the user prescribed function in cases
  of a linear initialization of L_0 E-M fields).

Particle initialization

• the number of particles in each cell of the physical domain EXPECT_EQ the user
  supplied number of particles per cell, on L0.

• On all levels : patch ghost particles are exact copies of particles in counterpart
  cells on neighbor patches of the same level.

• L0 has no levelGhostOldParticles, no levelGhostNewParticles and no levelGhostParticles

• L1 domain particles are equal to particles obtained from splitting L0 particles and that fall into L1 patches domain

• Splitting all particles of level L_i and keeping those in the level L_{i+1}
  (more refined) particle ghost region results in an particle array exactly equal
  to the levelGhostOldParticles and levelGhostParticles of level L_{i+1}

Moment initialization

in the following, we define as "complete node":

• all physical domain nodes

• all ghost node (patchGhost or levelGhost) that are supposed to receive a
  complete contribution from surrounding particles.

• On all levels, on each 'complete node', the density of each
  population equals the user supplied density up to a certain precision that
  varies as 1/sqrt(N) with N the number of particles
  per cell of the considered level. This is checked domain nodes.

• On all levels, on each 'complete node', the density of the ions
  EXPECT_DOUBLE_EQ the sum of the density
  over all populations. This is checked on domain nodes.

• on all levels, on each 'complete node', the flux of each population should be
  equal (up to a certain precision) to the product of the user supplied density
  and the bulk velocity of that population.
  user supplied bulk velocity components should all be different.

• on all levels, the bulk velocity of the ions should be the equal to the sum of
  all population fluxes divided by the ion density.

Goal

This feature will allow PHARE to refine the mesh at location satisfying a set of criteria based on the current state of the simulation.

we probably don't need to keep this https://github.com/PHAREHUB/PHARE/blob/master/src/core/utilities/function/function.h
because we already have those in :

https://github.com/PHAREHUB/PHARE/blob/master/src/initializer/data_provider.h

core depends on initializer (not the other way around) so probably we can keep those in data_provider and remove those in core?

For now it only loads particles but also needs to load electromagnetic fields.
Particles are loaded via a ParticleInitializer which comes out a factory which is given the particleInfo from each population. The particleInfo is the (sub)dict from the initialization.

Probably the same technique should be used for EM fields. EM fields Ctor would take a dict. containing its names, initialization procedures, etc.

splitting

In order to generalize the splitting to 2d & 3d, we need first to define the way to store all the weights & delta values. The 3 arrays below are associated to dimension = 1, 2 & 3. For each interpolation order (from 1 to 3 in 3 associated columns), w = m and d=n means we need mdifferent values of weights and n different values of deltas.

For exact splitting, we use capital letters : W & D. Empty squares, mean that this number of babies is not relevant, either because exact splitting would need less particles, or because we never calculated it with SplitPIC

dimension 1 :

dimension 2 :

dimension 3 :

• list all quantities, algorithms and schedules needed for coarsening and when
• how should MultiPhysicsIntegrator::standardLevelSynchronization use the Messengers for coarsening ?
• how can we use QuantityCommunicator for coarsening
• coarsening of the electric and magnetic field
• coarsening of the ion moments
• what is a good test of the coarsening done by the Messenger if the coarsening Operator is already test ?

built with

cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$PWD ..

fix - add to CMakeLists.txt

include_directories(
    ${CMAKE_CURRENT_SOURCE_DIR}/build/include
)

need to test periodicity on level 0 for:

• magnetic field components
• electric field components
• moments (ions and populations)
• particles (patch ghost == correct domain particles)
• no levelghost particles

for:

• all interp orders
• dimension 1 (2 and 3 can be prepared and fail for now)

User initialization of the root level finishes with calling fillRootGhosts.
This test intends to check that electromagnetic, moments and patchGhostParticles are well filled.

Definition

Fail fast/Fail early. 
Boundary testing and validation of python Simulation inputs.

Assumptions & Requirements

Certain sets of inputs are valid, others are not and should not validate.
Some inputs are related, and require cross-validation. i.e. max_nbr_levels has to be equal to or larger than the highest level set from "refined_boxes"

Needs

Assess the boundaries and validation appropriate for each simulation input value.
Execution test, with two sets of input Simulation objects with varying values
Extendable setup to add new valid/invalid versions as developments proceed.

Tests

Execute Simulation initialization based on various inputs
- Those that should pass
- Those that should fail

Questions:

Do we implement "proper_nesting_buffer" in tandem? It could impact the validation of the positions of refined_boxes I think.

case in point https://github.com/PHAREHUB/PHARE/blob/master/src/amr/data/field/field_data.h#L87

PHARE should be able to run and read the python initialization script expected in the PythonDataProvider.
Python initialization should be working with PHAREIN

We store our schedules in the messengers. Whenever SAMRAI removes a level, the associated schedules should be deleted as well.

Coarsening test WithRatioFrom2To10TestThat/ALinearFieldCoarsen1DO1.conserveLinearFunction/0 is broken with MPI.

• do we actually need that test if sub-components are tested?
• can we test the coarsening is OK for all quantities and all interp_order and all dimensions in a functional test?

Come back on that once #58 is done

seems the cmake file is using https://github.com/monwarez/SAMRAI

why not https://github.com/LLNL/SAMRAI ?

regrid is not tested yet and should be.
The difficulty here is to create a hierarchy in a configuration where regrid is needed so that isRegriddingis true in initializeLevelData and the routine is called.

Regrid affects all quantities they should all be tested.

moveIons is the function that takes electromagnetic fields, ions and move the ion distribution and moments.

The implementation can be inspired from the design project : https://hephaistos.lpp.polytechnique.fr/rhodecode/GIT_REPOSITORIES/LPP/phare/tests/particledesign/files/257aad76884cbc1fe142917bd82b1b5470002752/interface/solverppc.cpp#L134

currently we include both ./src and ./src/$MODULE

it's decided to just use ./src , and remove ./src/$MODULE per dependency

refinement_boxes": {"L0": {"B0": [(10,), (14,)]}}

doesn't work with

mpirun -n 2

it only starts to work when 14 is changed to 19

This test has been disabled for parallel execution.
The problem is that it assumes level1 only has 1 patch. In parallel, more patches are created and it is not easy to get the expected particles lying only in the levelGhostArea.

This test should be deleted once #62 is done

tests/amr/data/field/copy_pack/stream_pack/*

has some peculiar issues going on so the tests have been disabled in fix for #51

This ticket is to assess/fix/re-enable these tests.

Context

Level ghost particles contributing to the moments have two origins:

• those refined from the next coarser at its former time step : LevelGhostParticlesOld
• those refined from the next coarser at its current time step : LevelGhostParticlesNew

when calculating the moments on a given level at an arbitrary time, let alpha be the ratio of the current time since the last level synchronisation time (next coarser former time step) and the total duration of the next coarser time step, then the contribution of the level ghost particles is :

alpha * levelGhostParticlesNew + (1-alpha) * levelGhostParticlesOld.

alpha being 0 before the first time step of the sub cycle, and 1 at the end of the sub cycle

In ion_updater.h we have :

double alpha = 0.5;
interpolator_(std::begin(pop.levelGhostParticlesNew()),
                      std::end(pop.levelGhostParticlesNew()), pop.density(), pop.flux(), layout,
                      /*coef = */ alpha);


interpolator_(std::begin(pop.levelGhostParticlesOld()),
                       std::end(pop.levelGhostParticlesOld()), pop.density(), pop.flux(), layout,
                       /*coef = */ (1. - alpha));

This means that for now, alpha is hard-coded to be 0.5.

Objective

This issue aims at getting alpha defined correctly, i.e. to be the current time elapsed since last sync point. divided by the total duration of the next coarser time step.

hint:

this is done in hybrid_hybrid_messenger_strategy.h

        double timeInterpCoef_(double const beforePushTime, double const afterPushTime)
        {
            return (afterPushTime - beforePushTime)
                   / (afterPushCoarseTime_ - beforePushCoarseTime_);
        }

and used in the same file:

        virtual void fillIonMomentGhosts(IonsT& ions, SAMRAI::hier::PatchLevel& level,
                                         double const beforePushTime,
                                         double const afterPushTime) override
        {
            auto alpha = timeInterpCoef_(beforePushTime, afterPushTime);


            for (auto patch : level)
            {
                auto dataOnPatch = resourcesManager_->setOnPatch(*patch, ions);
                auto layout      = layoutFromPatch<GridLayoutT>(*patch);

                for (auto& pop : ions)
                {
                    // first thing to do is to project patchGhostParitcles moments
                    auto& patchGhosts = pop.patchGhostParticles();
                    auto& density     = pop.density();
                    auto& flux        = pop.flux();

                    interpolate_(std::begin(patchGhosts), std::end(patchGhosts), density, flux,
                                 layout);


                    // then grab levelGhostParticlesOld and levelGhostParticlesNew
                    // and project them with alpha and (1-alpha) coefs, respectively
                    auto& levelGhostOld = pop.levelGhostParticlesOld();
                    interpolate_(std::begin(levelGhostOld), std::end(levelGhostOld), density, flux,
                                 layout, 1. - alpha);

                    auto& levelGhostNew = pop.levelGhostParticlesNew();
                    interpolate_(std::begin(levelGhostNew), std::end(levelGhostNew), density, flux,
                                 layout, alpha);
                }
            }
        }

Regriding changes the hierarchy and the invalidate schedules.
The schedules listed in Refiners and Synchronizers must thus be updated after regriding. This is done in the resetHierarchyConfiguration() method of the StandardTagAndInitStrategy (in our case the MultiPhysicsIntegrator)

for now computeIonMoments_ is a private method of the HybridHybridMessengerStrategy. It is used in:

• fillRootGhosts()

it is not fillRootGhosts business to also compute ion moments, even if after fillRootGhosts() is called, we want moments OK on the grid.

• initLevel()
• regrid()

Later, computing ion moments from all particles will also be needed in the solver, in particular in moveIon(). Thus we need to move computeIonMoments_ outside the HybridHybridMessengerStrategy.

• HybridHybridMessengerStrat still needs to use it
• MultiPhysicsIntegrator cannot call computeIonMoments_() since it does not know what ion moments are. It only has messengers and models.

• Pharein should be a subproject for PHARE
• a python script ran by phare as it starts should :
  • take the python script building the Simulation object with all diagnostics and initial condition, etc.
  • add data from the Simulation object and use pyphare 'add(key, value)' to add the data to the PHAREDict.

configure GitHub and team city so they work together when PRs are created

PHARE needs to output data to disk.

• Write particles for each population
• Write density for ions and populations
• Write bulk velocity and flux for ions and populations
• Write Electromagnetic fields
• Datasets are written per patch per level
• both particle and field ghost nodes are written
• All levels specified by the user need to written at sync points.
• User can specify the dump times
• Diags are written in HDF5 using HighFive
• Dual/Primal value in each direction should be given for each field written or the Diagnostic file should specify the layout type (e.g. Yee) and python post-proc routines know what is where.
• write time of the dump as an HDF5 group (first after root)
• Diagnostics are specified by blocks in job.py and passed to pyphare in Simulation/Diagnostics/
• User can ask for either E or B or both in Electromag Diag
• User can ask for specific Fields or VecFields of total ions or specific Populations in FluidDiagnostics

https://github.com/PHAREHUB/PHARE/blob/master/tests/core/utilities/point/test_point.cpp

put global and per test coverage for PR builds

Refactor tests that use a BasicHierarchy or related class to use the parametrised hybrid hierarchy developed in #58
Define ranges for varying parameters related to these tests (interp_order, dimension, etc.)

Writing all particles to the disk is heavy and most of the time useless because rarely the user wants to see them all.

This issue is about adding to the pharein ParticleDiagnostics block, and to the cpp particle diagnostic, the possibility of selecting only specific particles, based on a criteria provided by the user.

the most general is the user providing a python function that gets contiguous particle vectors and returns a mask or a subset associated with the selected particles.

It would be convenient to get the python stack trace from pybind, when phare crashes due to errors in init.py.

This will be useful for setting up tests
We should be able to choose:

• interp_order (compile time parameter)
• dimension (compile time param)
• max number of levels (SAMRAI)
• min or max number of cells per patch (SAMRAI)
• and all parameters relevant for the hybrid simulation
• params should be configured from the python script
• user input refinement boxes for specific level

Use full hybrid hierarchy done in #58 to test all schedules in HybridMessenger work ok.
This includes schedules for initialisation and ghost filling, for fields and particles, for refining and synchronization.

This should fix #60, #59, #61

is this https://github.com/PHAREHUB/PHARE/blob/master/src/core/numerics/moments/moments.h#L33 correct ?
should we deposit levelGhostOld particles with domain and patch ghosts ?
levelGhostOld are defined at the first time step of the level subcycle (or last next coarser step). So it's fine to deposit these particles when initialising a level... but probably not when advancing the solution since during subcycling these particles are not at the current time and fillIonMomentGhost should rather be used (it uses both levelGhostNew and levelGhostOld)

The previous test cases which were changed to use the simulator and "job.py", only run job.py once.

This may not strictly be an issue as the variance is over the dimension and interp_order, and in these cases the explicit template parameters of the Simulator class are controlled from C++.

there should be a build conf (or several) on the CI that builds PHARE for versions of gcc from 7 to 9, and clang from 6 to 9. Use the already packaged versions of the compilers.

Field streaming test is disabled for parallel execution.
Rather than fixing it : Test streaming works OK in a functional test that uses all schedules PHARE is supposed to run. This should use all schedules used in the HybridMessenger.
Come back on that once #58 is done.

from some preliminary testing it seems the latest ICPC has issues with pybind11

to testing periodically

Not yet Implemented or Tested

Not yet implemented in 2D or 3D

1. In src/amr/wrapper/integrator.h, the integrator class only works at 1 dimension.

2. The refinement of particles in src/amr/data/particles/refine/split.h only works for few refinements, and not 3D.

3. I think in src/amr/data/field/coarsening/field_coarsen.cpp, only the 1D case is considered.

4. The electrons in src/core/data/electrons/electrons.h are only 1D.

5. There is no implementation of src/core/data/ions/ions.h at 2D & 3D.

6. Python file src/phare/phare_init.py is only 1D

Not yet tested in 2D or 3D

FieldData

• field variable test : needs to define variables living on and inside borders in 2 and 3D. Either transform the test in Typed_test to loop over dims, but loose the "parametric" aspect, or do one Test_P per dim... #411

• field time interpolate tests is only 1D. Generalizing for 2D and 3D is relatively straightforward and just consists in updating loops that fill and inspect the fields. (issue #365)

• AFieldLinearRefineIndexesAndWeights1D: are testing the startIndex and the Weights for refinement are working ok for dual and primal nodes, and for odd and even refinement ratios, all this in 1D. We do need to add these tests for 2D and 3D but only for refinement ratio 2 (thus only even). #412

• ALinearFieldRefine1DO1: This tests initializes Patch hierarchies where each quantity is initialized on level 0 with a linear function, that is refined for finer levels. The test consists in checking that the values on refined levels "double equal" the evaluation of the linear function at the level node coordinates, which is expected since the refinement operator is linear. This test should be generalized to 2D, 3D and also to all interporders. However we do not have to make these new versions working for refinement ratios 2 to 10, only 2 is enough. (issue #363)

• ALinearFieldCoarsen1DO1 This test uses a 2 levels basic hierarchy where each quantity is intialized on level 1 to a linear function and checks that after applying a coarsening schedule, values on the coarse level are correct. This test is only 1D. It is also failing in MPI probably because it assumes that user refinement boxes are the patches on the refined level... This test could be re made in python by initializing a hierarchy, advancing a time step and dumping diags. Advancing 1 coarse step would make a sync of the refined level and we could check coarse values where there are fine patches more easily from a python hierarchy. once done in python the cpp test can be deleted.

• AFieldCoarsenOperator: the test uses a python script to generate expected fine and coarsened values and checks the coarse values are as expected. The test is only 1D and needs to pass in 2D and 3D. The python script is ugly and could use some refactoring as well. 2D is now done, #413 for 3D

ParticleData

• AParticlesData1D (copy and stream) is only 1D and should be made 2D and 3D versions. (issue #364)

• levelOneCoarseBoundaries, areCorrectlyFilledByRefinedSchedule: this test checks that levelGhosts on the refined level of a 2 levels hierarchy are ok by splitting the whole L0 and filtering those which are in the ghost area box. Since there is no easy way to get only levelGhost boxes of a patch in C++ this only works when there is 1 patch on L1, or several but not touching. As a result the test only works in serial because parallel creates several Patch. This test should be deprecated by the #349 which implements the levelghost test (#286) see #414

• levelOneInterior, isCorrectlyFilledByRefinedSchedule: this tests that interior particles of refined levels are indeed exactly equal to the split of all particles of the next coarser restricted to the L1 domain box. The test is 1D and 2D. It could easily be 3D if an input_3d_ratio_2.txt was made. One could think re-doing this test as part of the initialization python tests since it is more flexible. see #415

Messengers

• HybridBasicHierarchyTest, fillsRefinedLevelFieldGhosts is 1D and 2D. Not sure we need to implement it in 3D given the test does not work with MPI : should be removed once #410 is done

Diagnostics

• Diagnostics are tested in the test Simulator2dTest and Simulator1dTest, the tests initializes the hierarchy, dump and read the file and compares it with the actual data in the hierarchy. This needs to be done in 3D as well. The tests are ran for all interp orders and all refined particle nbrs possible. see #416

• diagnostics dumps are tested from python. The test checks that a dump/advance()/dump() sequence indeed makes a file that contains both times (0 and first advance). All diags files are checked. The test runs for all interp orders but not all dimensions. Considering making it 2D and 3D is possible but the dimension should not matter for what's tested so it may be low priority. #275

DataWrangler

• The DataWrangler is tested for all interporders but only for 1D. Extention to 2D and 3D is necessary.

Splitting

• refined_particle_nbr test...?

Initializer

• test_initializer.cpp tests properties of the DataProviders, including that the PythonDataProvider can read a python script and provide a valid PHAREDict. This only is done for 1D and should be extended to 2D and 3D.

AMaxwellianParticleInitializer1D

This tests that a MaxwellianParticleInitializer loads the expected number of particles and that they are in the domain. It is only 1D. It could be generalized to 2D and 3D but I wonder if that is not already covered by the initialization tests.

14. The tests in tests/core/data/electrons/test_electrons.cpp are only 1D and with interporder = 1.

15. The test in tests/amr/models/test_main.cpp is only at 1D and interporder = 1.

16. The test in phare/pharein/examples/job.py is only 1D.

17. The test in tests/core/data/gridlayout/deriv.py is only 1D.

18. The test in tests/core/data/gridlayout/gridlayout_amr.cpp is only 1D.

19. The test in tests/core/data/maxwellian_particle_initializer/test_main.cpp is only 1D.

20. The test in tests/core/numerics/ampere/test_ampere.py is for 1, 2 & 3D, bu only for interporder=1

21. The test in tests/core/numerics/ion_updater/test_updater.cpp is only 1D

22. The test in tests/core/utilities/particle_selector/test_main.cpp is only 1D.

23. The test in tests/core/utilities/partitionner/test_main.cpp seems to be only 2D... not 1D neither 3D.

24. The tests in tests/initializer/job.py re only 1D.

25. The tests in tests/initializer/test_initializer.cpp are only 1D.

26. The tests in tests/simulator/py/init.py are only 1D.

27. The tests in tests/simulator/py/refinement_boxes.py are only 1D.

28. The tests in tests/simulator/per_test.h are only 1D.

Misc

1. The test in phare/pharein/simulation.py accept interpOrder, up to 4 (shouldn't it be 3 ?)
2. Max refinementFactor = 10 is hard coded in src/amr/data/field/coarsening/field_coarsen_operator.h

Most objects in phare are parametrised by compile time parameters.
The most basic of these parameters are interp_order and dimension. Most other parameters derive from those.
Although they are compile time, we need PHARE to let the user choose these parameters at runtime.

Using something like this :

https://hephaistos.lpp.polytechnique.fr/rhodecode/GIT_REPOSITORIES/LPP/phare/tests/solvercompiletime/files/251d737e1bb92ce919fc780199295087f80b8860/solvert.cpp

should make this possible.

Also see if this code can be simplified.