UWGeodynamics velocitySurface_2D Bug

Raising awareness of bugs within the velocitySurface_2D function in UWGeo that @julesghub and I have identified last week in both UWV2.13 and V2.15.

My usage:

 Model.surfaceProcesses = GEO.surfaceProcesses.velocitySurface_2D( 
     airIndex=air.index, sedimentIndex=sediment.index,
     sedimentationRate= 1. *u.millimeter / u.year, erosionRate= 1. * u.millimeter / u.year,
     surfaceElevation=0.*u.kilometer,
     surfaceArray = coords_surface)

In V2.13, the function operates correctly depositing sediments at the defined rate, but upon restart sediments instantly replace all the air below the defined surface elevation (0 km in my case) in one timestep. This behaviour means we cannot restart models while using this sedimentation function.

The animation shows the model swarms evolution with a model restart 6 seconds in. Basement is coloured green, sediments are yellow, and air is grey.

https://github.com/underworldcode/underworld2/assets/50604575/aad34db0-640f-46f1-9ac4-d53f1d7ae23c

In V2.15, models will not restart with the sedimentation function. Attempting to restart a model with this function enabled returns this error:

Error

``` --------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In[7], line 10 2 Model.surfaceProcesses = GEO.surfaceProcesses.velocitySurface_2D( 3 airIndex=air.index, sedimentIndex=sediment.index, 4 sedimentationRate= 0.5 *u.millimeter / u.year, erosionRate= 0.5 * u.millimeter / u.year, 5 surfaceElevation=0.*u.kilometer, 6 surfaceArray = coords_surface) 8 Model.init_model() ---> 10 Model.run_for(110000 * u.years, checkpoint_interval=20000. * u.year, restartStep=-1, restartDir="Vel_Surf_O") 11 # %% File /opt/venv/lib/python3.10/site-packages/underworld/UWGeodynamics/_model.py:1657, in Model.run_for(self, duration, checkpoint_interval, nstep, checkpoint_times, restart_checkpoint, dt, restartStep, restartDir, output_units) 1652 ndduration = self._ndtime + nd(duration) if duration else None 1654 output_dt_units = _get_output_units( 1655 output_units, checkpoint_interval, duration) -> 1657 checkpointer = _CheckpointFunction( 1658 self, duration, checkpoint_interval, 1659 checkpoint_times, restart_checkpoint, output_dt_units) 1661 if not nstep: 1662 nstep = self.stepDone File /opt/venv/lib/python3.10/site-packages/underworld/UWGeodynamics/_model.py:2396, in _CheckpointFunction.__init__(self, Model, duration, checkpoint_interval, checkpoint_times, restart_checkpoint, output_units) 2393 self.outputDir = Model.outputDir 2395 if checkpoint_interval or checkpoint_times: -> 2396 self.checkpoint_all() File /opt/venv/lib/python3.10/site-packages/underworld/UWGeodynamics/_model.py:2477, in _CheckpointFunction.checkpoint_all(self, checkpointID, variables, tracers, time, outputDir) 2475 self.checkpoint_fields(variables, checkpointID, time, outputDir) 2476 self.checkpoint_swarm(variables, checkpointID, time, outputDir) -> 2477 self.checkpoint_tracers(tracers, checkpointID, time, outputDir) 2478 comm.Barrier() File /opt/venv/lib/python3.10/site-packages/underworld/UWGeodynamics/_model.py:2673, in _CheckpointFunction.checkpoint_tracers(self, tracers, checkpointID, time, outputDir) 2671 if Model.passive_tracers: 2672 for (dump, item) in Model.passive_tracers.items(): -> 2673 item.save(outputDir, checkpointID, time) 2675 comm.Barrier() File /opt/venv/lib/python3.10/site-packages/underworld/UWGeodynamics/_utils.py:193, in PassiveTracers.save(self, outputDir, checkpointID, time) 191 obj.data[...] = field["value"].evaluate(self) 192 else: --> 193 obj.data[...] = field["value"].evaluate(self.data) 195 handle = obj.save('%s.h5' % file_prefix, units=field["units"]) 197 if rank == 0: 198 # Add attribute to xdmf file File /opt/venv/lib/python3.10/site-packages/underworld/function/_function.py:938, in Function.evaluate(self, inputData, inputType) 936 if (not (inputData.base is None)) or inputData.flags['F_CONTIGUOUS']: 937 inputData = inputData.copy() --> 938 return _cfn.Query(self._fncself).query(_cfn.NumpyInput(inputData,inputType)) 939 else: 940 # try convert and recurse 941 return self.evaluate( self._evaluate_data_convert_to_ndarray(inputData), inputType ) RuntimeError: Issue utilising function of class 'SwarmVariable' constructed at: 0- 0:/tmp/ipykernel_36/2600013166.py:2 Model.surfaceProcesses = GEO.surfaceProcesses.velocitySurface_2D( Error message: Function input dimensionality (2) does not appear to match swarm variable dimensionality (1952542047). ```

underworldcode / underworld2

UWGeodynamics velocitySurface_2D Bug #675