Creating a new code

Hi, I am a new user of gpusph and I’ve some queries about how to modify the given codes in problems to better suite my requirement. I am trying to modify the code by changing the rheology to herschel_bulkley but where do I provide the yield strength? For changing the input parameters such as height of channel, time step, and dimensions, should I change it directly in .cu file ? Moreover I would like to exploit the output to obtain the velocity characteristics of the flow; y_v curve but I can’t understand the output generated.

Hi @Ayan_Ganguly

Concerning the parameters of the simulation, yes, you need to change the .cu file where you put all the parameters you like. Then compile GPUSPH and execute !

For the post processing, all this can be done in paraview without any problems. Just extract the fields you want from plot data.

regarding the first the point ( rheology to herschel_bulkley), I am sorry I dont experience that. Surely, you will have a reply from my colleagues that experience better this point.


Hello @Ayan_Ganguly and thanks for your interest in GPUSPH.

To set the yield strength for non-Newtonian fluids, you can use set_yield_strength(fluid_idx, ys) where fluid_idx is the fluid index (that gets returned by add_fluid; if it’s a single-fluid simulaiton, you can just hard-code the value 0) and ys is the value of the yield strength. For the exponent, you can use set_visc_power_law(fluid_idx, n) with n being the power to which the strain rate should be raised.

The output on console is only used as a quick reminder of the progress of the simulation, the actual data is stored in a directory under tests/, usually named after the problem name and the date/time when the simulation was launched. You can open the saved data found in that directory (.vtp files) with ParaView, and plot any of the fields you care about.

Please do keep in mind that not all non-Newtonian rheological models have been validated yet (Papanastasiou and Bingham should work, for the others we haven’t run validation tests, since the very high viscosity found in the non-yielding region benefit significantly from the use of a semi-implicit integrator, which has only been recently completed).

I’m also working on a Python script (that still requires ParaView) to extract specific field data, but it’s not in the repository yet. I hope to find the time to merge it soon.

Thanks @giuseppe.bilotta, I’m an undergrad student and I would like a bit more clarification please. Where are the functions set_yield_strength(fluid_idx, ys) and set_visc_power_law(fluid_idx, n) defined because I can’t get the auto completion prompt for either function in the file. Also where is the round_up function defined?
I get that the add_fluid() is used to pass the density of the fluid as described in the setup pdf but in OpenChannel test file the function takes an input argument of 2650.0f; why so?
And Sir finally, whenever I define the still water level around 3X mm, I get an warning as follows:

which furthermore generates a bizarre peak velocity. How to overcome this?

Hello @Ayan_Ganguly,

the set_yield_strength and set_visc_power_law functions are defined in ProblemCore.h. I’m not sure why your editor doesn’t suggest it for completion. The round_up function is defined in utils.h.

The reason why OpenChannel sets the fluid density to 2650 kg/m³is because that’s the density of lava, which inspired the test case. For water one would use1000 instead, of course.

Setting the still-water level affects the hydrostatic filling of the particles and the maximum fall velocity, and thus also the sound speed (if not set manually). This can affect the stability of the simulation when it’s too low.

Sir, I’ve been trying the simulation for non-newtonian rheology with the Herschel-Bulkley model but every time there is a FATAL message during runtime after the first iteration requesting to quit simulation and it ends as:

For this I’ve edited the framework model as :
I’ve a couple of questions regarding this as in;
When does this message get prompted and what are the modifications required?
Thank you

Hello @Ayan_Ganguly,

the fatal error occurs because the time-step required to satisfy the CFL constratints has become too small. This can happen for two reasons:

  1. a possible reason is that the simulation is unstable and this leads to unphysical results that ultimately cause this kind of errors. A common cause for this is when the user sets the time-step manually at the beginning: even if this is only used during the first iteration (if dynamic time-stepping is enabled), it may still be problematic. One can in general trust the automatic computation of the time-step, so commenting any set_timestep in the code, or setting it to a low but not extremely low value (e.g. 1.0e-6) may help. Among the messages shown when the problem starts execution there is information about the auto-computed time-step requirements, and whether or not the initial user choice fits with them or not.

  2. another possible reason, which can manifest with extremely low Reynold numbers or with non-Newtonian fluids that form a plug (and thus effectively infinite viscosity) is that the time-step requirement for the viscous term are so stringent that for a stable simulation the time-step simply does need to be that low. The only practical option in this case (at least until support for the semi-implicit integration scheme gets integrated into next which should hopefully happen soon-ish) is to try and run the simulation at a lower resolution: since the time-step constraint from the viscous term scales with the square of the resolution, you may not even need to reduce the resolution (set_deltap) too much.

I hope this helps identify and resolve the issue in your test case (without additional information about the setup, it’s impossible for me to tell exactly what is going wrong). Do let us know if this is the case, or if you need further assistance.


Giuseppe Bilotta