Test: 2D base case

[bartvanwesten]

We need to implement tests for the following scenario:

Description

Testing the transport-related functionalities step-wise for a 5x5 2D grid without supply-limitations. The test includes starting, running and finishing the model, checking the generation of output files and the actual values of generated output. The main aim is to use as many as defaults as possible, to test reproducability / backwards-compatibility.

Code Reference

[ Permalink to code for which the test is required ]

Test Cases

Case:

• A 5x5 2D grid (x and y from 0-4), non-rotated grid with a bed level at z=0, without bed-update and default non-erodible layer (?)

• Wind series for 36 timesteps, with increasing wind direction (0-360 deg) in Cartesian direction

• dt = 100, tstart = 0 and tstop = 3600 (=36*100)

• Bagnold formulation, 1 sediment fraction with a grainsize of 250 mu

• No masks, vegetation, hydrodynamics or meteo

• Only processes are wind, transport and threshold

• Only threshold is grainsize

• All boundaries are constant

• Outputtime = dt

• Output parameters are zb, uw, uws, uwn, ustar, ustars, ustarn, uth, Cu (and maybe Ct, q, qs, qn)

• Preconditions:

• ...

• Test Steps:

• Start, run and finish the model and find, read and check the output

• Test Data:

• Input files, see above

• Expected Result:

• Model runs and finishes

• Model generates output-file

• Checks related to output file:

  • Does the output-file exist?
  • Amount of parameters in output file?

• Checks related to the output parameters

  • dimension of parameter (10, 5, 5, ...) ?
  • value of parameter (check with hard coded)

• Postcondition: [ conditions that need to achieve when the test case was successfully executed ]

[Sierd]

Test should ideally include a check for continuity. Which means for any domain:

• total erosion = outbound flux + sediment in transport or
• cum_sum(zb) = cum_sum(qs_total(at boundaries)) + sum(Ct(at t_end))

This gets complicated when having more grain sizes than 1.

[niketagrawal]

Implementation of this test is being done here. Changes are work in progress.

Suggestions for testcases to include:

• Check for garbage values such as 'NaN' and '-1' in the netCDF file.
• Check the input file for a valid and realistic configuration. For example, wind speed of 100 m/sec is not realistic. This could be the first test case to run in the for loop of the integration test. It is desired that if such values are encountered in the input file, a warning must be issued on the console, instead of an error. This behavior could be implemented in the source code.
• Run the model for a particular case by turning off/on a combination of certain process flags. This could be achieved by using one file as the base file and creating versions of it with desired process flags switched Off/ON which would be used as input for the subsequent run.

[niketagrawal]

Changing the solver and time step did not influence the mass deficit. Possible changes needed in the code @Sierd

[niketagrawal]

To do: Test cases to check basic functionality of the model:

• Compare entire netCDF files. Content is mostly deterministic.
• Number, name and shape of parameters as defined in aeolis.txt should be the same in the netCDF file and they should not be empty.

[Sierd]

@niketagrawal, I have drafted a test for continuity. Could you try and run this and see if this works. The idea is that frac is a very small number at the end of the simulation (the plot shows frac for every timestep).

This seems to work for me, can you double check?

import netCDF4
import numpy as np
import matplotlib.pyplot as plt

fname = r'.\examples\1D\case1_small_waves\aeolis.nc'

with netCDF4.Dataset(fname, 'r') as ds:

        # dimensions
        dt = np.diff(ds.variables['time'][:])[0]
        dx = np.diff(ds.variables['x'][0,:])[0]
        zb = ds.variables['zb'][:,:,:]
        qs = ds.variables['qs']        
        Ct = np.sum(ds.variables['Ct'][:,:,:,:], axis=-1)

        steps = zb.shape[0]-1
        dV=np.zeros(steps)
        frac=np.zeros(steps)

        for i in range(0,steps):
            # bed level change
            dz = zb[i,:,:] - zb[0,:,:]
            V1 = dz * dx #* dy

            # volume loss over downwind border
            V2 = qs[:i,0,-1,0]  * dt / (2650. * .6)

            # volume in saltation
            V3 = Ct[i,:,:] * dx / (2650. * .6)

            # metrics
            E = np.abs(np.minimum(0., V1).sum())
            D = np.abs(np.maximum(0., V1).sum())
            dV[i] = np.abs(V1.sum() + V2.sum() + V3.sum())
            frac = dV/E

plt.plot(frac)
plt.show()