Coupling AeoLiS and OpenFoam

[Sierd]

In this thread I would like to collect experiences related to the coupling of AeoLiS and OpenFOAM. I get many emails about this and I will summarize a bit of the experience so far.

Coupling AeoLiS and OpenFOAM would make sense because:

1. AeoLiS' fast sediment transport calculations could extent the results of OpenFOAM simulations towards sediment transport and morphological changes.
2. OpenFOAM can provide a much more elaborate calculation of winds over topography and the consequences for shear stress near the bed than the somewhat rudimentary (but very vast) interpretations in AeoLiS. AeoLiS' sediment transport calculations may benefit from a more detailed calculation of wind shear.

Dr. Pourteimouri and co-workers have successfully set up a coupling of OpenFOAM and AeoLiS and applied it to study different configurations of beach infrastructure.
https://www.sciencedirect.com/science/article/pii/S0169555X23002118

Also, Dr. Cecil and co-workers are working on a coupling
https://doi.org/10.5194/egusphere-2024-855

Both examples above used the "hotstart" functionality but other ways of interfacing may be possible.

The interface between OpenFOAM and AeoLiS may be facilitated by:

1. The BMI functionality of AeoLiS is made to exchange information between different models. It would make sense to exchange ‘ustar’ from OpenFOAM to AeoLiS and to exchange ‘zb’ from AeoLiS to OpenFOAM. This exchange will be subject to timeintervals that depend on calculations times etc. Technically this could be done. The BMI functionality can easily be used using a callback function.
   https://aeolis.readthedocs.io/en/latest/user/implementation.html#basic-model-interface-bmi
   After looking at this in more detail it appears that AeoLiS currently overwrites the ‘ustar’ after the callback so there needs to be a change to facilitate the use of BMI to do this.
2. The hotstart functionality was initially built to hotstart models that take a long time after they crashed for some reason during their calculations. Hotstart files can be used to transfer OpenFOAM results to AeoLiS. The documentation is not very mature. The Hotstart functionality can be found here.
   

   aeolis-python/aeolis/model.py

   Line 3119 in 6afe8ee

   def load_hotstartfiles(self) -> None:

3. ...
   (I would be curious to learn about alternatives.)

I would like to explore the best way of coupling AeoLiS and OpenFoam. Please add your thoughts on this.

[saeb-faraji-gargari]

I tested importing the shear stress calculated by OpenFOAM into AeoLiS using a BMI script (as taus), but it does not seem to work as expected. The shear stress is initialized with a zero array instead of using the imported OpenFOAM shear stress. The AeoLiS setup and the OpenFOAM shear stress files are attached.

The OpenFOAM shear stress files are named wallShearStress_0_grid.grd, wallShearStress_1_grid.grd, and wallShearStress_2_grid.grd for the x (s), y (n), and z directions, respectively.

I would appreciate any guidance.
aeolis_setup_saeb.zip

[Sierd]

Hi Saeb,

Thank you for your input. I ran your model, it seems to run on my computer as well. I did not check output yet. I understand that in your callback.py function you try to set 'taus'.

I understand that the callback for updating tau is not working. I think this is because the callback runs before the update function. This happens in this loop.

aeolis-python/aeolis/model.py

Lines 2825 to 2835 in 6afe8ee

	# start model loop
	self.t0 = time.time()
	self.output_write()
	while self.t <= self.p['tstop']:
	if callback is not None:
	callback(self)
	self.update()
	self.output_write()
	self.print_progress()

In the update function the taus gets overwritten by the AeoLiS model calculations. I have confirmed in my tests that this is indeed the case.

Now I think about it, I realize that using the callback function as it is programmed now, only works when adapting the cumulative variables (like 'zb' and model parameters like 'V_ver'). Instantaneous variables get overwritten in the model. This is an issue that we can solve from a technical perspective by reading the callback on another place.

However before doing that, I think it is interesting to look at this from a more physical perspective. You are trying to calculate sediment transports. You are using the 'bagnold' formulation. In the code that can be found here.

aeolis-python/aeolis/transport.py

Lines 351 to 357 in 6afe8ee

	if p['method_transport'].lower() == 'bagnold':
	s['Cu'][ix] = np.maximum(0., p['Cb'] * rhoa / g * (ustar[ix] - uth[ix])**3 / u[ix])
	s['Cuf'][ix] = np.maximum(0., p['Cb'] * rhoa / g * (ustar[ix] - uthf[ix])**3 / u[ix])
	s['Cu0'][ix] = np.maximum(0., p['Cb'] * rhoa / g * (ustar0[ix] - uth0[ix])**3 / u[ix])

This is where an equilibrium concentration in a cell is calculated by estimating the transport (using a standard Bagnold formulation based on ustar) divided by the average velocity of the grains. The latter is a very complicated procedure based on Duran's work.

I wonder if the details provided by OpenFOAM can help us a bit. Can OpenFoam give us ustar and an estimate for grainspeed (equal to some percentage of the near-bed velocity or something like that). This could simplify some of the assumptions related to grainspeeds. I think it would be interesting to make an 'OpenFOAM' version of the sediment transport calculation in AeoLiS that is tailored to the coupling with OpenFOAM. This way we can experiment with this coupling a bit in a more isolated way.

Let me know what you think. I would be happy to help coding such a process.

[saeb-faraji-gargari]

Hi Sierd,

Thank you for your response.

OpenFOAM computes the shear stress on the ground, and the corresponding u_star can be calculated as u_star = sqrt(tau / rho), where tau is the shear stress and rho is the fluid density.

[Sierd]

Excellent Saeb,

We can use ustar to estimate sediment transports using Bagnold.

Does it also give some output on near bed velocity? Like 20cm above the bed? My idea is to estimate sediment velocities based on their extremes between 0 and the near bed velocity as a maximum to start with.

I think we can use ustar to calculate sediment transport (using Bagnold) and convert that to concentrations by dividing by the sediment velocity. This division by the sediment velocity is needed for the advection solver which is the basis of AeoLiS. I feel that this way we allow for a relatively light and seamless connection between OpenFOAM and AeoLiS without overlapping complexities.

If you agree we can start doing the technical work. I created a branch where we can work in:
https://github.com/openearth/aeolis-python/tree/OpenFoam_patch

I made already a keyword here:

aeolis-python/aeolis/constants.py

Line 191 in 47ad11b

'OpenFOAM' : False, # NEW! [T/F] Boolean to indicate coupling with OpenFOAM.

(False is the default, we set it True in our input file.)

The keyword is not functional yet. We may want to deactivate the wind/shear solvers in AeoLiS and activate the coupling with OpenFOAM. We can use input files generated by OpenFOAM to start with.
My first idea is to use BMI and make the callback function operational for this purpose. We can also hard-code some stuff first to test a bit.

What do you think?

[saeb-faraji-gargari]

Hi Sierd,

Thank you for sharing your idea.

Extracting the velocity above the ground surface (for instance, a few centimeters) is possible, although it requires processing and interpolation from the OpenFOAM results. It is important to note that OpenFOAM computes the air velocity, so the sediment-air velocity (u_sed) needs to be calculated based on this.

[saeb-faraji-gargari]

One additional point:

If we use the modified Bagnold equation (Nickling and Neuman, 2009 see Equation 1 in https://www.sciencedirect.com/science/article/pii/S0169555X23002118, the ground wall shear stress can be used to calculate the sediment flux (q), which is then used to compute c_eq (see Equation 5 in the mentioned paper). In this way, I think we do not need to extract u_star a few centimeters above the bed, although we still need to extract the wind velocity for u_sed.

[saeb-faraji-gargari]

Hi Sierd,

I would like to open an issue to work on further coupling between OpenFOAM and Aeolis.

see [issue #229] [https://github.com//issues/229]