Gridap & Richards equation

Specific Domains Modelling & Simulations
1

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
2

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
3

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
4

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
5

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!

7

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 const keyword 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_everywhere function 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