Numerical questions about spatially varying parameters for ParallelStencil.jl

Hi developers for ParallelStencil.jl,

Hope you stepped into the new year with ease. I have experienced with your julia software in a while, thanks to your clear documentation and workflow!

I have a couple of questions about ParallelStencil.jl/HydroMech2D.jl at main · omlins/ParallelStencil.jl · GitHub
while I am using the software:

  1. How do we allow the background porosity distribution to be spatially varying? Right now the background porosity is set to be homogenous as $$\phi_0$$ and only accepts scalar (since other physical parameters also depend on it, e.g. time step). Same for density, etc.
  2. When simulating with larger background porosity (e.g. 0.2) and larger perturbation (e.g. 0.4), the behavior does not look like chimneys but more like “blobs”. Is this behavior expected? Also, do I have to change the physical time stepping and reduction time stepping accordingly? Is there anything I should take care of in terms of numerical stability of the algorithm?
  3. Also, out of curiosity, will the adjoint-based inversion solver be available in this ParallelStencil.jl package in the future?

Thanks a lot for your help!

Kind Regards,


Hi Francis, thank you for getting in touch with @samo and myself regarding the HydroMech2D.jl featured example from the ParallelStencil.jl package. Also thanks for the feedback about documentation.

Regarding your questions:

  1. Background porosity is initialised in the Phi array, as defined here; Then an elliptical perturbation is added. If you want to perturb the background (now constant) porosity, you could define a different initialisation (e.g. adding random noise) replacing line 147. For density, both fluid and solid density are defined as scalar, the total (Phi-averaged) density being an array. If you want to implement spatially variable e.g. solid density (times gravity), you could define a new array in line 127 such as Rosg = @ones(nx ,ny ). Then you could add your perturbation to Rosg, pass it as Data.Array::Rosg to the compute_params_∇!(...) kernel and change following line as such: @all(Rog) = ρfg*@all(Phi) + @all(Rosg)*(1.0-@all(Phi)) - ρgBG

  2. Larger background porosity values trigger larger effective permeability values. The effective permeability ultimately controls the compaction length, which is the characteristic length scale of the two-phase system (e.g. Räss et al. 2019). So yes, this behaviour is expected. Pseudo-transient time-stepper Vsc, Ptsc, Pfsc and physical time step dt_red scaling factors may need to be adapted accordingly.

  3. Yes, the 2D and 3D adjoint-based inversion solvers (based on Reuber et al. 2020) are current work in progress and should be published as well as available as ParallelStencil.jl miniapp in the near future.

Hope this helps, and thanks for your interest !


Hi Ludovic,

Thanks for your quick reply!

  1. Background porosity also affects solid shear viscosity in line 114. It would therefore affect the calculation of time stepping. More importantly, qDy would be affected as well. Could you please give me some instructions on how to fix them to accept the heterogenous background model?

  2. Could you also give me some brief instruction on adapting time stepping? e.g. if I double the background porosity, how should I change these time stepping? Also, if we change the porosity to be heterogenous, permeability will also be heterogenous. How do we calculate the compaction length then?

  3. That sounds great!!!

Thanks again for your support!


Hi Francis,

Could you please give me some instructions on how to fix them [μs, qDy] to accept the heterogenous background model?

For the shear viscosity, you could to following:

  • comment out line 114;
  • define a “host” array for μs: Mus = zeros(nx ,ny ) after line 145;
  • initialise Mus with your newly implemented spatially variable background porosity array Phi0, or with another function of Phi, e.g. Mus .= ηC0*Phi0/η2μs;
  • copy the host array to the device array Mus = Data.Array(Mus).

Then one will have to modify all the places where the scalar (Data.Number) μs was used, replacing it by the Data.Array array Mus. E.g. kernel definitions will have to take Data.Array::Mus as argument now.

Since the shear viscosity is now an array, you will have to change the definition of the pseudo-time steps. You can to comment out lines 120-121, initialise 3 array:

dτVx = @zeros(nx-1,ny-2)
dτVy = @zeros(nx-2,ny-1)
dτPt = @zeros(nx  ,ny  )

and the following kernel to compute them:

@parallel function compute_dτ!(...)
    @all(dτVx) = min_dxy2/@av_xi(Mus)/(1.0+β_n)/4.1/Vsc
    @all(dτVy) = min_dxy2/@av_yi(Mus)/(1.0+β_n)/4.1/Vsc
    @all(dτPt) = 4.1*@all(Mus)*(1.0+β_n)/max(nx,ny)/Ptsc

taking care to pass the correct arguments to the new kernel (here ...). See the 2D Stokes miniapp for reference.

The y-direction Darcy flux qDy calculation you point to is the vertical boundary condition definition of the inflow and outflow of the model. The definition is chosen to satisfy an effective pressure Pt-Pf=0 at the boundaries to avoid compaction/decompaction at the boundaries. If you want to have non-scalar values there, you could perform vectorised calculation on host arrays or create a kernel to precompute a vector for the lower and upper boundary. This may actually not be really needed.

Could you also give me some brief instruction on adapting time stepping […] how should I change these time stepping

There is not much to do there. You can play with dt_red, the controls on the physical time step reduction, to see if you better resolve the waves depending on the configuration.

Also, if we change the porosity to be heterogenous, permeability will also be heterogenous. How do we calculate the compaction length then?

Porosity is heterogenous. What you change is the background porosity, the “reference value” of the variable porosity. Permeability is and will further be heterogenous as it depends on Phi and not on on ϕ0. To compute the compaction length, you can choose which reference permeability you want to consider. I would take the far field, background (or reference) one.

I hope this helps. If you need further support, feel free to start a dev git repo; I could then provide more direct support for your particular needs and project.


1 Like

Thanks for your instructions, Ludovic! I’ve forked your repo and created a spatially varying porosity example at here. I overloaded a bunch of parallel functions as your instruction. Now it works without any error. However, there is a problem that some of the intermediate variables (Rx, Ry, Vx, Vy etc) get to NaN in only a couple of iterations, which finally blows out the simulation.

In this example, I was trying to set up a shale layer in the middle region, and sandstone layers on the upper and lower sides. The parameters (bulk viscosity, permeability) follow your table 2 in your nature paper. Since I choose the shale layer to be the reference one, the compaction length of the model is quite small so lx, ly are chosen to be very large to make a model size in the order of kilometers.

Could you please help me see what’s wrong in the simulation? Thanks for your support again!