# Gridap & Richards equation

**URL:** <https://discourse.julialang.org/t/gridap-richards-equation/123961>\
**Category:** Modelling & Simulations\
**Tags:** gridap\
**Created:** [December 18, 2024, 2:57pm UTC](https://discourse.julialang.org/t/gridap-richards-equation/123961 "2024-12-18T14:57:01Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Fabio\_Zottele](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabio_zottele/32/214136_2.png) [@Fabio\_Zottele](https://discourse.julialang.org/u/Fabio_Zottele)\
**Post date:** [December 18, 2024, 2:57pm UTC](https://discourse.julialang.org/t/gridap-richards-equation/123961/1 "2024-12-18T14:57:01Z")

</div>

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](https://en.wikipedia.org/wiki/Richards_equation)). 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…

---

<div class="post-metadata">

**Author:** ![Paulo\_Jabardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paulo_jabardo/32/3196_2.png) [@Paulo\_Jabardo](https://discourse.julialang.org/u/Paulo_Jabardo)\
**Post date:** [December 18, 2024, 5:35pm UTC](https://discourse.julialang.org/t/gridap-richards-equation/123961/2 "2024-12-18T17:35:18Z")

</div>

I believe it is suitable. Checkout the tutorial examples

[https://gridap.github.io/Tutorials/stable/](https://gridap.github.io/Tutorials/stable/)

In particular checkout the 4th example and the 7th.

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [December 19, 2024, 9:09pm UTC](https://discourse.julialang.org/t/gridap-richards-equation/123961/3 "2024-12-19T21:09:27Z")

</div>

> [@Fabio\_Zottele](#):
>
> I would like to know if Gridap is suitable for solving “highly non linear problem” (the Richards equation, [Richards equation - Wikipedia](https://en.wikipedia.org/wiki/Richards_equation)).

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:

> **[A posteriori error estimates and adaptivity in numerical approximation of PDEs...](https://dumas.ccsd.cnrs.fr/THESES-SU/tel-04631602v1)**
>
> This thesis concerns a posteriori error analysis and adaptive algorithms to approximately solve nonlinear partial differential equations (PDEs). We consider PDEs of both elliptic and degenerate parabolic type. We also study adaptivity in floating...

---

<div class="post-metadata">

**Author:** ![aerappa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aerappa/32/212281_2.png) [@aerappa](https://discourse.julialang.org/u/aerappa)\
**Post date:** [December 20, 2024, 12:39am UTC](https://discourse.julialang.org/t/gridap-richards-equation/123961/4 "2024-12-20T00:39:05Z")

</div>

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](https://gridap.github.io/Tutorials/dev/pages/t018_transient_nonlinear/). 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.

---

<div class="post-metadata">

**Author:** ![Fabio\_Zottele](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabio_zottele/32/214136_2.png) [@Fabio\_Zottele](https://discourse.julialang.org/u/Fabio_Zottele)\
**Post date:** [December 21, 2024, 8:43pm UTC](https://discourse.julialang.org/t/gridap-richards-equation/123961/5 "2024-12-21T20:43:15Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![aerappa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aerappa/32/212281_2.png) [@aerappa](https://discourse.julialang.org/u/aerappa)\
**Post date:** [December 24, 2024, 10:58pm UTC](https://discourse.julialang.org/t/gridap-richards-equation/123961/7 "2024-12-24T22:58:30Z")

</div>

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:

```julia
θ(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

```julia
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

```julia
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!