# Numerical questions about spatially varying parameters for ParallelStencil.jl
**URL:** https://discourse.julialang.org/t/numerical-questions-about-spatially-varying-parameters-for-parallelstencil-jl/53558
**Category:** Modelling & Simulations
**Tags:** package
**Created:** [January 18, 2021, 8:06pm UTC](https://discourse.julialang.org/t/numerical-questions-about-spatially-varying-parameters-for-parallelstencil-jl/53558 "2021-01-18T20:06:38Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@ziyiyin97](https://discourse.julialang.org/u/ziyiyin97)
#### Post date: [January 18, 2021, 8:06pm UTC](https://discourse.julialang.org/t/numerical-questions-about-spatially-varying-parameters-for-parallelstencil-jl/53558/1 "2021-01-18T20:06:38Z")
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](https://github.com/omlins/ParallelStencil.jl/blob/main/miniapps/HydroMech2D.jl)
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,
Francis
---
### Author:  [@luraess](https://discourse.julialang.org/u/luraess)
#### Post date: [January 19, 2021, 10:08am UTC](https://discourse.julialang.org/t/numerical-questions-about-spatially-varying-parameters-for-parallelstencil-jl/53558/2 "2021-01-19T10:08:25Z")
> [@ziyiyin97](#):
>
> ParallelStencil.jl/HydroMech2D.jl at main · omlins/ParallelStencil.jl · GitHub
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](https://github.com/omlins/ParallelStencil.jl/blob/8d6b32efe9a682713960d4724a7dbb556ed115bf/miniapps/HydroMech2D.jl#L147); Then an elliptical perturbation is [added](https://github.com/omlins/ParallelStencil.jl/blob/8d6b32efe9a682713960d4724a7dbb556ed115bf/miniapps/HydroMech2D.jl#L150). If you want to perturb the background (now constant) porosity, you could define a different initialisation (e.g. adding random noise) replacing [line 147](https://github.com/omlins/ParallelStencil.jl/blob/8d6b32efe9a682713960d4724a7dbb556ed115bf/miniapps/HydroMech2D.jl#L147). 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](https://github.com/omlins/ParallelStencil.jl/blob/8d6b32efe9a682713960d4724a7dbb556ed115bf/miniapps/HydroMech2D.jl#L25) 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](https://doi.org/10.1093/gji/ggz239)). 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](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwisq5D036fuAhVP-6QKHR9mDH8QwqsBMAB6BAgEEAM&url=https%3A%2F%2Fwww.youtube.com%2Fwatch%3Fv%3D1t1AKnnGRqA&usg=AOvVaw3jTK0YnixBqUXnBasUJIpd) (based on [Reuber et al. 2020](https://doi.org/10.1016/j.jcp.2020.109797)) 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 !
Ludovic
---
### Author:  [@ziyiyin97](https://discourse.julialang.org/u/ziyiyin97)
#### Post date: [January 19, 2021, 3:24pm UTC](https://discourse.julialang.org/t/numerical-questions-about-spatially-varying-parameters-for-parallelstencil-jl/53558/3 "2021-01-19T15:24:55Z")
Hi Ludovic,
Thanks for your quick reply!
1. Background porosity also affects solid shear viscosity in [line 114](https://github.com/omlins/ParallelStencil.jl/blob/8d6b32efe9a682713960d4724a7dbb556ed115bf/miniapps/HydroMech2D.jl#L114). It would therefore affect the calculation of time stepping. More importantly, [qDy](https://github.com/omlins/ParallelStencil.jl/blob/8d6b32efe9a682713960d4724a7dbb556ed115bf/miniapps/HydroMech2D.jl#L154) 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!
Francis
---
### Author:  [@luraess](https://discourse.julialang.org/u/luraess)
#### Post date: [January 19, 2021, 9:45pm UTC](https://discourse.julialang.org/t/numerical-questions-about-spatially-varying-parameters-for-parallelstencil-jl/53558/4 "2021-01-19T21:45:34Z")
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:
```Julia
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:
```julia
@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
return
end
```
taking care to pass the correct arguments to the new kernel (here `...`). See the [2D Stokes miniapp](https://github.com/omlins/ParallelStencil.jl/blob/8d6b32efe9a682713960d4724a7dbb556ed115bf/miniapps/Stokes2D.jl#L11) for reference.
The y-direction Darcy flux `qDy` calculation you [point to](https://github.com/omlins/ParallelStencil.jl/blob/8d6b32efe9a682713960d4724a7dbb556ed115bf/miniapps/HydroMech2D.jl#L154) 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.
Ludovic
---
### Author:  [@ziyiyin97](https://discourse.julialang.org/u/ziyiyin97)
#### Post date: [January 20, 2021, 4:44am UTC](https://discourse.julialang.org/t/numerical-questions-about-spatially-varying-parameters-for-parallelstencil-jl/53558/5 "2021-01-20T04:44:42Z")
Thanks for your instructions, Ludovic! I’ve forked your repo and created a spatially varying porosity example at [here](https://github.com/ziyiyin97/ParallelStencil.jl/blob/main/miniapps/HydroMech2D_vary.jl). 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!
Francis