Add TimeStepper base class plus ForwardBackward and RK4 time steppers

[mwarusz]

This PR adds a base class for all Omega time steppers and two concrete implementations of simple explicit time steppers:

• forward-backward scheme
• fourth order Runge Kutta scheme

A test for the above is also added. This PR depends on #109 and includes it with some modifications. Since this PR is missing comments and documentation and is based on a draft PR I am marking it also as a draft.

@hyungyukang

Checklist

• Documentation:
  • User's Guide has been updated
  • Developer's Guide has been updated
  • Documentation has been built locally and changes look as expected
• Testing
  • CTest unit tests for new features have been added per the approved design.
  • Unit tests have passed. Please provide a relevant CDash build entry for verification.

[hyungyukang]

Added the RK2 time stepper and opened a PR to mwarusz#1.

[mwarusz]

I merged @hyungyukang's RK2 PR into this branch.

[mark-petersen]

Passes all tests on perlmutter with CPU and GPU:

details:

######### perlmutter CPU
CODE_DIR=opr

cd /global/homes/m/mpeterse/repos/e3sm/${CODE_DIR}
git submodule update --init --recursive externals/ekat externals/scorpio cime
cd components/omega/

module load cmake
RUNDIR=240819_omega_pr_c
mkdir ${PSCRATCH}/runs/$RUNDIR
cd !$

rm -rf build
mkdir build
cd build
export PARMETIS_ROOT=/global/cfs/cdirs/e3sm/software/polaris/pm-cpu/spack/dev_polaris_0_4_0_nvidia_mpich/var/spack/environments/dev_polaris_0_4_0_nvidia_mpich/.spack-env/view

cmake \
   -DOMEGA_CIME_COMPILER=nvidia \
   -DOMEGA_BUILD_TYPE=Release \
   -DOMEGA_CIME_MACHINE=pm-cpu \
   -DOMEGA_PARMETIS_ROOT=${PARMETIS_ROOT}\
   -DOMEGA_BUILD_TEST=ON \
   -Wno-dev \
   -S /global/homes/m/mpeterse/repos/E3SM/${CODE_DIR}/components/omega -B .
./omega_build.sh

# linking:
cd test
ln -isf /global/homes/m/mpeterse/meshes/omega/O*nc .


# run test:
salloc --nodes 1 --qos interactive --time 01:00:00 --constraint cpu --account=m4572 # or e3sm
RUNDIR=240819_omega_pr_c
cd ${PSCRATCH}/runs/${RUNDIR}/build

./omega_ctest.sh

# output in Testing/Temporary/LastTest.log

######### perlmutter GPU
CODE_DIR=opr
RUNDIR=240819_omega_pr_gpu_b
mkdir ${PSCRATCH}/runs/$RUNDIR
cd !$

rm -rf build
mkdir build
cd build
module load cmake

export PARMETIS_ROOT=/global/cfs/cdirs/e3sm/software/polaris/pm-gpu/spack/dev_polaris_0_4_0_nvidiagpu_mpich/var/spack/environments/dev_polaris_0_4_0_nvidiagpu_mpich/.spack-env/view
cmake \
   -DOMEGA_CIME_COMPILER=nvidiagpu \
   -DOMEGA_BUILD_TYPE=Release \
   -DOMEGA_CIME_MACHINE=pm-gpu \
   -DOMEGA_PARMETIS_ROOT=${PARMETIS_ROOT}\
   -DOMEGA_BUILD_TEST=ON \
   -Wno-dev \
   -S /global/homes/m/mpeterse/repos/E3SM/${CODE_DIR}/components/omega -B .
./omega_build.sh

# linking:
cd test
ln -isf /global/homes/m/mpeterse/meshes/omega/O*nc .

# run test:
salloc --nodes 4 --qos interactive --time 01:00:00 --constraint gpu --tasks-per-node=2 --gpus-per-task 1 --account=m4572_g # or e3sm_g
RUNDIR=240819_omega_pr_gpu_b
cd ${PSCRATCH}/runs/${RUNDIR}/build

./omega_ctest.sh

If you copy these instructions, change paths to your own.

[mark-petersen]

@mwarusz thank you for the additional comments within the code. That is helpful!

[hyungyukang]

To fully test this PR, I quickly built a test version of the shallow water model using this PR. It was relatively straightforward to integrate this PR with the previous Omega-0 test model by invoking doStep from this PR.

I have also added the struct for the manufactured solutions (ManufacturedThicknessTendency{}, ManufacturedVelocityTendency{}), with initial condition files generated by Polaris. Additionally, I turned on tendency computations manually (bool Enabled = false => true) in TendencyTerms.h, which is fine for now as Omega-0 has a capability reading a YAML config file.

In the test model, the time-stepping component is structured as follows (as a rough example):

   Config Options;

   auto *TestAuxState =
       AuxiliaryState::create("TestAuxState", Mesh, NVertLevels);
   auto *TestTendencies =
       Tendencies::create("TestTendencies", Mesh, NVertLevels, &Options,
                          ManufacturedThicknessTendency{}, ManufacturedVelocityTendency{});

   auto *TestTimeStepper = TimeStepper::create(
       "TestTimeStepper", TimeStepperType::RungeKutta2, TestTendencies, TestAuxState, Mesh, DefHalo);

   const auto *Stepper = TimeStepper::get("TestTimeStepper");

   for (int step = 0; step < nsteps; ++step) {
       // Compute current time
       CurrentTime = step * dt;
       Stepper->doStep(State, CurrentTime, dt);
       sw_check_status(printInterval, CurrentTime, dt, Comm, Env, Mesh, State);
   }

With the different time steppers implemented in this PR (ForwardBackward, RungeKutta2, RungeKutta4), I obtained the following convergence rates across various horizontal resolutions:

Based on the results above, it can be concluded that the time steppers implemented and the other changes introduced in this PR are functioning correctly. @mwarusz , thank you very much for your efforts on this PR!

I have additional minor comments, so I will add them to the codes soon.

[mark-petersen]

Thank you @mwarusz for you dedicated effort on this PR. Approving based on visual inspection, my testing on perlmutter, and @hyungyukang's convergence plot and testing by adding RK2.

[sbrus89]

This is looking great. I just have a couple minor comments so far.

[philipwjones]

Since I'll be heading out on vacation, I will go ahead and approve now based on a code read-through and testing by others. The code is already starting to get somewhat more complex than I'd like so the added comments and a good DevGuide section will be important.

[hyungyukang]

Based on @mark-petersen 's ctests on Perlmutter, my own full functional testing with the manufactured solution and visual inspection by others, I approve this PR. The other changes look good as well. @mwarusz , thank you very much for your hard work on this PR.

If you have any other concerns about the time level changes, please disregard my suggestions.

[hyungyukang]

I assume tendencies for the manufactured solution will be computed here?

[sbrus89]

@mwarusz, @hyungyukang - Do we want to add the manufactured solution tendency terms to this PR, or add them in another?

[mwarusz]

Let's do it in another PR.

[hyungyukang]

I was wondering if it would be easy to import some constants defined in a YAML file here. For example, the wavenumber of the manufactured solution is defined in the YAML file and we need to include it here. I think it should be easy, but I thought it would be good to check.

[mwarusz]

It should be straightforward. Just make the constructor of this class take in Config.

[brian-oneill]

I'm a little confused by the Time being included as an argument of the doStep function here, and in the various compute tendencies functions in the Tendency class. The tendency calculations are not generally dependent on the time directly, and the only thing I'm seeing it used for in the PR is for computing the exact solution in the unit test.

[brian-oneill]

Also, since the TimeStep is kind of core to the purpose of the TimeStepper, maybe makes sense to just store it as a member variable of the class

[mwarusz]

The manufactured solution has forcings that depend explicitly on time. I agree about the TimeStep. Initially I thought that making it an argument will make it easier to support varying it during a simulation, but now I think that, if we want to support that, having a changeTimeStep method would be better.

[sbrus89]

For down the road, tidal potential forcing also depends on time.

[brian-oneill]

Okay, I was thinking of the TimeStepper as being more or less equivalent to the various time integrator subroutines in mpas-ocean, which only directly take the time step as an argument, but I guess in mpas they have the overall simulation time buried in modules for sim-time-dependent forcings (like ocn_vel_tidal_potential_tend). When we end up adding forcings we could could have the simulation time stored within class objects representing those forcings, but this is all kind of arbitrary and can be worked out as we go along.

I do think that as opposed to a Real representing Time, it should probably be a TimeInstant object, which represents an instant in time (obviously) on a particular Calendar object, which it contains a pointer to. There is a NoCalendar type, in which the TimeInstant just becomes the elapsed time in seconds, which can be used for unit tests and test problems like the manufactured solution.

[mwarusz]

I made the TimeStep a member, and changed Time and TimeStep to use types TimeInstant and TimeInterval, respectively.

[brian-oneill]

Ran the unit tests successfully. Aside from maybe some quibbles on structure, looks good to me.

[sbrus89]

I was running into issues with a local merge into develop, but with the rebase, this now passes CTests on Frontier CPU and GPU. Thanks for all your work on this @mwarusz!

[mwarusz]

I re-ran the tests and everything passes on Chrysalis, Frontier CPU and GPU, and Perlmutter CPU and GPU.

[mark-petersen]

Thanks! I will review documentation today.

[sbrus89]

One comment I have is that it may make more sense for each Timestepper to set the value for NTimeLevels as opposed to it being an independent option in the YAML input file. However, I'm not sure the best way to do that given that the State is initialized outside of the Timestepper.

[mwarusz]

One comment I have is that it may make more sense for each Timestepper to set the value for NTimeLevels as opposed to it being an independent option in the YAML input file. However, I'm not sure the best way to do that given that the State is initialized outside of the Timestepper.

@sbrus89 I see two options:

• The TimeStepper is initialized before the State and the default state gets this value from the default time stepper.
• We create a function that maps the time stepper config option to the TimeStepperType enum and a mapping from the TimeStepperType enum to NTimeLevels. The the state init doesn't need the time stepper to be initialized, it only needs to know which timestepper the user chose.

[sbrus89]

@mwarusz, I think the first option is the most straightforward. Perhaps we could have an optional argument to the State constructor that requires you either propvide an explicit NTimeLevels value, or a TimeStepper instance? That way we don't have a hard dependency on the TimeStepper being initialized.

[mwarusz]

@sbrus89 I implemented the first option in 44268f8

[sbrus89]

@sbrus89 I implemented the first option in 44268f8

That looks good to me, thanks @mwarusz

[hyungyukang]

After several updates in this PR, I have tested this branch again by implementing a shallow water model with a time stepping loop from this PR. I conducted two tests:

1. Global steady state nonlinear geostrophic flow on a sphere (Williamson et al., 1992)

• L2 thick error norm remains steady under 1e-4 (dt=100 s, dx = 240 km), which is similar to Ringler et al. (2010).
• Total energy is conserved well through the simulation.

2. nonlinear manufactured solution test on a Cartesian plane (Bishnu et al., 2024).

• The model demonstrates second-order accuracy with varying horizontal resolutions using both RK2 and RK4 time-stepping methods.

Based on the above results, I can say the most recent version of this branch is working properly. Thanks @mwarusz again for your work!

[mark-petersen]

The developers doc is good! I just had the one question.