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
return
end
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.
Ludovic