pwolfram / mpas-model Goto Github PK

@pwolfram, I think the dimensions of verticalTreatment should be changed ti just nParticles from (nParticles, Time). Would I be correct in thinking a particle should never transform its mode? I.e., it is only dependent on particle integer.

This becomes an issue when post-processing. In the following example (with trimmed down output):

ds = xr.open_dataset('mpaso.hist.am.lagrPartTrack.0034-01-01_00.00.00.nc')
ds.attrs = ''
print(ds)

<xarray.Dataset>
Dimensions:            (Time: 15, nParticles: 1037200)
Dimensions without coordinates: Time, nParticles
Data variables:
    indexToParticleID  (nParticles) int32 ...
    particleDIC        (Time, nParticles) float64 ...
    verticalTreatment  (Time, nParticles) int32 ...
    xParticle          (Time, nParticles) float64 ...
    xtime              (Time) |S64 ...
    yParticle          (Time, nParticles) float64 ...
    zLevelParticle     (Time, nParticles) float64 ...
    zParticle          (Time, nParticles) float64 ...

If I want to select for just surface floats to slim down the size:

ds = ds.where(ds.verticalTreatment==1, drop=True)
print(ds)

<xarray.Dataset>
Dimensions:            (Time: 15, nParticles: 64825)
Dimensions without coordinates: Time, nParticles
Data variables:
    indexToParticleID  (nParticles, Time) float64 9.724e+05 ... 1.037e+06
    particleDIC        (Time, nParticles) float64 1.923e+03 ... 2.189e+03
    verticalTreatment  (Time, nParticles) float64 1.0 1.0 1.0 ... 1.0 1.0 1.0
    xParticle          (Time, nParticles) float64 1.78e+06 ... -3.444e+06
    xtime              (Time, nParticles) object b'0034-01-02_00:00:00                                             ' ... b'0034-01-30_00:00:00                                             '
    yParticle          (Time, nParticles) float64 -4.414e+06 ... 8.133e+05
    zLevelParticle     (Time, nParticles) float64 -2.637 -1.952 ... -2.345
    zParticle          (Time, nParticles) float64 4.235e+06 ... 5.298e+06

note that indexToParticleID and xtime get converted to different data types and attempt to broadcast to the (Time, nParticles) format of verticalTreatment. If I select just one time slice of verticalTreatment, the where function proceeds as expected:

ds = ds.where(ds.verticalTreatment.isel(Time=0)==1, drop=True)
print(ds)

<xarray.Dataset>
Dimensions:            (Time: 15, nParticles: 64825)
Dimensions without coordinates: Time, nParticles
Data variables:
    indexToParticleID  (nParticles) float64 9.724e+05 9.724e+05 ... 1.037e+06
    particleDIC        (Time, nParticles) float64 1.923e+03 ... 2.189e+03
    verticalTreatment  (Time, nParticles) float64 1.0 1.0 1.0 ... 1.0 1.0 1.0
    xParticle          (Time, nParticles) float64 1.78e+06 ... -3.444e+06
    xtime              (Time, nParticles) object b'0034-01-02_00:00:00                                             ' ... b'0034-01-30_00:00:00                                             '
    yParticle          (Time, nParticles) float64 -4.414e+06 ... 8.133e+05
    zLevelParticle     (Time, nParticles) float64 -2.637 -1.952 ... -2.345
    zParticle          (Time, nParticles) float64 4.235e+06 ... 5.298e+06

E3SM namelist / stream options need updated following merge of new LIGHT capabilities.

Free floating particles require specialized boundary conditions for the top and the bottom, especially due to partial bottom cells.

@pwolfram, I think it would be great to have an option to append the instantaneous velocity fields to lagrPartTrack output. I.e., the velocities that are interpolated to the particle position to advect them. In isopycnal particles, I imagine this would just be one value. In passive, maybe this is u,v,w or just magnitude.

Context: I'm reading through analysis methods for computing total transport from Lagrangian particles and Drake et al. 2018 (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017GL076045) mention weighting transport trajectories by instantaneous velocity multiplied by the area of the cell the particle is at (methods, final paragraph).

Under new developments with free floating and surface particles, ensure that outputed indexLevels are updated to correspond to each particle type.

@pwolfram ,

Currently pyAMG coarsens the initial mesh in 3 dimensions, whereas it should probably only coarsen in 2D and then apply this filter to all vertical layers. I.e., we should only use pyAMG to coarsen the horizontal so that the vertical seeding option remains valid in the vertical.

For instance, here's the linear seeding case for an arbitrary grid cell in the 60to30km grid:

And after applying pyAMG:

The odd thing is that it still seeds the vertical column with 20 particles, but they don't at all follow the seeding strategy (e.g. linear or log).

cc @pwolfram @maltrud

The goal here is to give a few options for seeding in the vertical, in addition to just seeding linearly to the bottom. I'll work on this and add updates here.

Ensure that recent code changes are not breaking LIGHT under high-resolution cases.

Particle output should be reported in latitude and longitude coordinates to facilitate global analysis.

Additional regression suite tests should be added for new particle behaviors to make sure correct behavior is maintained under the possibility of future refractoring.

A particle mixed layer entrainment parameterization should be implemented to handle mixed layer dyanmics, e.g., transfers of particles to and from mixed layer.

The current use of pyAMG in make_particle_file.py to coarsen particle seeding in the horizontal is unfeasible for high resolution meshes. From memory, coarsening by 1 level in the 60to30km case takes ~45min.

I tried to coarsen 1 level in the 30to10km case (i.e. ./make_particle_file.py -i input.nc -g graph.9600 -p 9600 -o particles.nc --spatialfilter SouthernOceanXYZ -t 'passive','surface' -d 1) and the script is still running ~20 hours later. Is this an inherent limitation of pyAMG?

See

cc @mukundraj

Following LIGHT improvements, the framework helper function mpas_get_vertical_id
should also be updated.

It is stored on branch https://github.com/pwolfram/MPAS-Model/tree/framework/cleanup_mpas_get_vertical_id

For default streams file in testing_and_setup/compass/utility_scripts/LIGHTparticles/streams.ocean,

And add under the "lagrPartTrackOutput" stream:

<var name="currentCellGlobalID"/>
<var name="lonParticle"/>
<var name="latParticle"/>

To make these standard output for particles

cc @bradyrx

config_AM_lagrPartTrack_timeIntegration=4 needs to be added as a default option in E3SM, as noted in #9.

Scalar interpolation, e.g., for temperature and salinity, should also include horizontal interpolation (e.g., using the vertex-based barycentric interpolation) in order to better resolve fronts. The current process of only vertically-interpolating is only sufficient at very large scales.

Supercycling particle integration implemented in MPAS-Dev/MPAS#1484 should be verified for global runs.

Per @pwolfram , super-cycling should already be available for particles.

What is the most straightforward way to set (and thus vary) the super-cycling time length?

I checked Registry_lagrangian_particle_tracking.xml and didn't find anything related to it. However, I did find the super-cycling section in mpas_ocn_lagrangian_particle_tracking.F.

The global Eulerian cell ID corresponding to the particle location should be output to facilitate easy analysis for postprocessing.

In global runs, RK4 should be the default option for accuracy since time stepping with super cycling will not be using the Eulerian timestep dt. Larger time steps, e.g., 4-6 hrs, will benefit from use of RK4 at the cost of doubling each Lagrangian diagnostic step by a factor of 2 with respect to computation and memory.

@pwolfram, per our email discussion.

It would be helpful if you could test https://github.com/bradyrx/automate_mpas_simulation on clusters that aren't grizzly or wolf.

The script has been run for 240km, 60to30, and 30to10 for physics and BGC. Perhaps it would be best to include as a utility in the branch you will be merging to MPAS (i.e. this one)?

After you give it a spin, we should talk about crucial features to be useful in its first implementation. I'm thinking I should add:

• Branch runs for ocean/sea-ice

• RK schemes and supercycling for particles

• Flexibility on compsets (currently only supports g-cases)

• Edit particle output for sub-daily

• Edit highFrequency output for new variables and frequency