Hi all. I am rather new to julia and totally new to Gridap. Before starting a new time and energy consuming project I would like to know if Gridap is suitable for solving “highly non linear problem” (the Richards equation, Richards equation - Wikipedia). I would solve this equation first in 1D and then in 2D but having multi-material (different soil horizon with different hydraulic properties). For now I am just writing down the 1D-one material problem (weak forms, jacobians etc) on a paper sheet and before starting to code I would like some indications by you…
2 Likes
I believe it is suitable. Checkout the tutorial examples
https://gridap.github.io/Tutorials/stable/
In particular checkout the 4th example and the 7th.
2 Likes
Yes it most definitely is suitable: @aerappa spent part of his PhD thesis doing this. See for example his thesis, especially chapter 3 and the papers referenced therein:
3 Likes
Hi @Paulo_Jabardo , as @ffevotte mentioned I was able to solve some benchmark tests of the Richards equation using Gridap.jl. However, I basically wrote the nonlinear solver and time integration myself because I needed certain intermediate quantities between nonlinear steps to compute a posteriori error estimators. In addition, there was no high level API at the time in Gridap.jl for nonlinear transient problems.
However, it appears that such an interface now exists as per Tutorial 18. I think this would probably be a good place to start, but if you have some specific questions don’t hesitate to send me a message.
3 Likes
Thanks @aerappa @Paulo_Jabardo @ffevotte : I am studying Tutorial 4th, 7th and 18th and I made some advance in understanding how Gridap works. I am now translating the weak form and the jacobians in julia (i am new to julia’s grammar too, so I hope in coding everything correctly). However I still have hard time in setting the (very simple) boundary and initial conditions. I’ll try to figure out myself how to do, but winter holidays are coming and I’ll just take a break from the Richards equation!
Dear @aerappa , @Paulo_Jabardo & @ffevotte :,
I post here my work on solving the Richards’ equation for a 1D homogeneous soil profile using Gridap. The domain is defined as follows:
\Omega = [0, -1] \quad \text{(m)}
where ( z = 0 ) corresponds to the soil surface, and ( z = -1 ) represents the bottom of the soil profile. The domain is discretized into 50 mesh elements. I am trying to set the following boundary and initial conditions:
-
Boundary Condition: A Dirichlet condition (constant head for all the simulation) is applied at the soil surface:
h(z = 0) = 0\quad \forall t \geq 0 -
Initial Condition: The head is constant throughout the domain at ( t = 0 ):
h(z < 0) = -10000 \, \text{(Pa)}, \quad t = 0
I have written the weak forms for the Richards’ equation, the Jacobians for both the diffusion and capacity terms, and the necessary soil hydraulic property functions based on the van Genuchten model. However, I am struggling to correctly implement these initial and boundary conditions using Gridap.
Here is my code so far:
using Gridap
using GridapGmsh
θ_s = 0.45 # Saturated volumetric water content
θ_r = 0.05 # Residual volumetric water content
α = 0.1 # van Genuchten parameter (1/m)
n = 1.5 # van Genuchten parameter
m = 1 - 1/n # Relation between m and n
K_s = 1e-4 # Saturated hydraulic conductivity (m/s)
λ = 0.5 # Empirical exponent for K(h)
g = 9.81 # Gravitational constant (m/s^2)
# Funzioni per il modello di van Genuchten
function Se(h)
h >= 0 ? 1.0 : (1 + (α * abs(h))^n)^(-m)
end
function θ(h)
θ_r + Se(h) * (θ_s - θ_r)
end
function K(h)
Se_h = Se(h)
K_s * Se_h^λ * (1 - (1 - Se_h^(1/m))^m)^2
end
function C(h)
h >= 0 ? 0.0 : α * n * (1 - m) * (α * abs(h))^(n - 1) * (θ_s - θ_r) * Se(h)^(1 + m)
end
function ∂Se(h)
if h >= 0
return 0.0
else
sign_h = h < 0 ? -1 : 1
factor = 1 + (α * abs(h))^n
return -α * n * m * (α * abs(h))^(n - 1) * factor^(-m - 1) * sign_h
end
end
function ∂K(h)
Se_h = Se(h)
∂Se_h = ∂Se(h)
if Se_h == 1.0
return 0.0
end
term1 = λ * Se_h^(λ - 1) * (1 - (1 - Se_h^(1/m))^m)^2
term2 = 2 * Se_h^λ * (1 - (1 - Se_h^(1/m))^m) * m * (1 - Se_h^(1/m))^(m - 1) * (1/m) * Se_h^(1/m - 1)
return K_s * (term1 + term2) * ∂Se_h
end
function ∂C(h)
if h >= 0
return 0.0
else
Se_h = Se(h)
∂Se_h = ∂Se(h)
term1 = (n - 1) * α^(n - 1) * abs(h)^(n - 2) * Se_h^(1 + m)
term2 = α^(n - 1) * abs(h)^(n - 1) * (1 + m) * Se_h^m * ∂Se_h
return α * n * (1 - m) * (θ_s - θ_r) * (term1 + term2)
end
end
function h_van_genuchten(θ)
Se = (θ - θ_r) / (θ_s - θ_r)
h = -1 / α * ((1 / Se)^(1/m) - 1)^(1/n)
return h
end
model = GmshDiscreteModel("./MESHES/homogeneus_soil_1D.msh")
order = 1
reffe = ReferenceFE(lagrangian, Float64, order)
V = TestFESpace(model, reffe, conformity=:H1, dirichlet_tags=["Top"])
h0_top = h_van_genuchten(θ_s)
h0_domain = h_van_genuchten(θ_r+0.01)
U = TrialFESpace(V, h_top)
diffusion_term(t, u, v) = ∫(K(u) * ∇(u + z) ⋅ ∇(v))dΩ
capacity_term(t, u, v) = ∫(C(u) * ∂t(u) * v)dΩ
source_term(t, q, v) = ∫(q * v)dΩ
q = 0
residuals(t, u, v) = diffusion_term(t, u, v) + capacity_term(t, u, v) - source_term(t, q, v)
jacobian_diffusion(u, v, du) = ∫((∂K(u) * ∇(u + z) ⋅ ∇(v) * du) + (K(u) * ∇(du) ⋅ ∇(v)))dΩ
jacobian_capacity(u, v, du, dtu) = ∫((∂C(u) * dtu * v * du) + (C(u) * v * dtu))dΩ
Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)
op = TransientFEOperator(residuals, (jacobian_diffusion, jacobian_capacity), U, V)
linear_solver = LUSolver()
Δt = 30
θ_method= 0.5
solver = ThetaMethod(linear_solver, Δt, θ_method)
tableau = :SDIRK_2_2
solver_rk = RungeKutta(linear_solver, linear_solver, Δt, tableau)
t0, tF = 0.0, 1.0*60*60*24
uh0 = interpolate_everywhere(h_top, U)
uh = solve(solver_rk, op, t0, tF, uh0)
I am seeking guidance on how to:
- Set the boundary condition h(z=0)=0 \quad\forall t≥0 . This shoul be yet implemented.
- Initialize h(z<0)=−10000 Pa for t=0. I tried to use the
interpolate_everywherefunction with no succes.
I would greatly appreciate your suggestions or examples on how to apply these conditions correctly in Gridap.
Thank you in advance for your help!
Hi @Fabio_Zottele ,
First of all, a few minor “best practices” for Julia:
-
Try to avoid global variables, but in your case since they are constants you can use the
constkeyword for them. This is for performance reasons but also for debugging it can be helpful -
If you have a one line function you can write it like this:
θ(h) = θ_r + Se(h) * (θ_s - θ_r)
Now for your questions:
Set the boundary condition h(z=0)=0,\quad \forall t\ge 0. This shoul be yet implemented.
Essential boundary conditions (Dirichlet in your case) are imposed using the trial space. You are trying to do this here
U = TrialFESpace(V, h0_top)
However, your problem is that h0_top is not a function, but rather the value of the function h_van_genuchten evaluated at the constant θ_s. Instead, you should just do something like
h0_top(x) = h_van_genuchten(x[1])
because Gridap expects an nD vector x (even though you’re just in 1D for now)
Initialize h(z<0)=−10000Pa for t=0. I tried to use the
interpolate_everywherefunction with no succes.
You are likely experiencing the same problem here: you are interpolating a constant function equal to h_van_genuchten(θ_s) for all x.
Hope this helps!
2 Likes