# Help with solving transient equation using GetFEM (GWFL approach)

**URL:** <https://getfem.discourse.group/t/help-with-solving-transient-equation-using-getfem-gwfl-approach/70>\
**Category:** General\
**Created:** [March 24, 2025, 8:24pm UTC](https://getfem.discourse.group/t/help-with-solving-transient-equation-using-getfem-gwfl-approach/70 "2025-03-24T20:24:14Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Djoumessi](https://avatars.discourse-cdn.com/v4/letter/d/da6949/32.png) [@Djoumessi](https://getfem.discourse.group/u/Djoumessi)\
**Post date:** [March 24, 2025, 8:24pm UTC](https://getfem.discourse.group/t/help-with-solving-transient-equation-using-getfem-gwfl-approach/70/1 "2025-03-24T20:24:14Z")

</div>

Dear GetFEM community,

I am new to GetFEM and currently in the process of learning how to use it effectively. I would like to solve the two-variable FitzHugh-Nagumo equation using GetFEM, specifically via the **Generalized Weak Formulation Language (GWFL)**. I have studied the wave equation example, which relies heavily on bricks, but I want to implement my problem using GWFL instead.

The FitzHugh-Nagumo system is given by:

\begin{aligned} \frac{\partial u}{\partial t} &= D\_u \Delta u + u(1 - u)(u - a) - v, \\ \frac{\partial v}{\partial t} &= \epsilon \left( u - \gamma\_1 v + \beta\_1 \right) \end{aligned}

Below is the code I wrote to solve this system using GWFL.

```auto
import numpy as np
import getfem as gf
import os

# Paramtres du probleme
L = 1.0 # size of the cube
NX = 20 # mesh szie
NZ = 3 # mesh szie
dt = 0.1 # time step
Tmax = 2.0 # Temps final

# FitzHugh-Nagumo parametera
Du = 0.01 # Coefficient de diffusion pour u
epsilon = 0.1 # Paramètre epsilon
gamma1 = 0.5 # Paramètre gamma
beta1 = 0.1 # Paramètre beta
a = 0.1 # Seuil d'activation

# mesh
m = gf.Mesh('cartesian', np.linspace(0, L, NX),
                     np.linspace(0, L, NX),
                     np.linspace(0, L/4, NZ))
m.optimize_structure() # Convertit en tétraèdres

# FE spaces
mf = gf.MeshFem(m, 1)
mf.set_classical_fem(1) # Élément P1
mfrhs = gf.MeshFem(m,1) # Right-hand side
mfrhs.set_classical_fem(1)

# integration
mim = gf.MeshIm(m, 4)    

'''
# initial excitation
pts = np.array(mf.basic_dof_nodes())
print("Type de pts:", type(pts))
print("Forme de pts:", pts.shape)

if pts.shape[0] < 3:
    raise ValueError("Le tableau pts n'a pas 3 lignes. Vérifiez le maillage ou basic_dof_nodes().")

# Extraction des coordonnées x, y, z
x, y, z = pts[0, :], pts[1, :], pts[2, :]

# gaussienne avec NumPy
U = np.exp(-100 * ((x - 0.5)**2 + (y - 0.5)**2 + (z - 0.5)**2))
V = np.zeros_like(U)'''

U = mf.eval('-1*((x-0.5)*(x-0.5)*(y-0.5)*(y-0.5)*(z-0.125)*(z-0.125))') # Impulsion initiale au centre
V = np.zeros(U.shape)

# Model
md = gf.Model('real')

# Variables
md.add_fem_variable('u', mf)
md.add_fem_variable('v', mf)
md.add_fem_data("Previous_u", mf)
md.add_fem_data("Previous_v", mf)

# Initialisation 
#md.add_initialized_fem_data('Previous_u', mf, U)
#md.add_initialized_fem_data('Previous_v', mf, V)

md.set_variable('Previous_u', U)
md.set_variable('Previous_v', V)

#Newmark scheme
'''beta = 0.25
gamma = 0.5
md.add_Newmark_scheme('u', beta, gamma)
md.add_Newmark_scheme('v', beta, gamma)
md.set_time_step(dt)'''

#diffusion terme
md.add_initialized_data('Du', [Du])
md.add_initialized_data("dt", dt)
##md.add_nonlinear_term(mim, 'Du*(Grad_u.Grad_Test_u)')
#md.add_nonlinear_term(mim, 'Du*Grad(u).Grad(Test_u)')
#md.add_generic_elliptic_brick(mim, 'u', 'Du')

# nonlinear part
md.add_initialized_data('a', [a])
md.add_initialized_data('gamma1', [gamma1])
md.add_initialized_data('beta1', [beta1])
md.add_initialized_data('epsilon', [epsilon])

md.add_nonlinear_term(mim, '(u-Previous_u)*Test_u -dt*Du*Grad(u).Grad(Test_u) - dt*Test_u.(u*(1 - u)*(u - a) - v)')
md.add_linear_term(mim, '(v-Previous_v)*Test_v - dt*Test_v . epsilon * (u - gamma1*v + beta1)')

# dirichlet but no flux is the one used
#md.add_Dirichlet_condition_with_multipliers(mim, 'u', mf, 1)
#md.add_Dirichlet_condition_with_multipliers(mim, 'v', mf, 1)

# results folder
if not os.path.exists('Resultf'):
    os.makedirs('Resultf')

# time loop
t = 0
n = 0
while t < Tmax:
    print(f"Temps {t:.2f}")
    
    # Solve
    md.solve()
    
    # extraction
    U = md.variable('u')
    V = md.variable('v')

    # Export in VTK
    mf.export_to_vtk(f'Resultf/u_{n:04d}.vtk', U)
    mf.export_to_vtk(f'Resultf/v_{n:04d}.vtk', V)
  
    # Update variable
    md.set_variable('Previous_u', U)
    md.set_variable('Previous_v', V)
    
    t += dt
    n += 1
    #md.shift_variables_for_time_integration()

print("Simulation terminée. Les résultats sont exportés en VTK.")

```

However, when running it, I encounter the following error messages:

```auto
Temps 0.00
mbly RHS
Trace 2 in getfem_models.cc, line 2653: Global generic assembly RHS
Trace 2 in getfem_models.cc, line 2653: Global generic assembly RHS
Trace 2 in getfem_models.cc, line 2655: Global generic assembly tangent term

```

I suspect the issue might be related to the RHS and Tangent. However, despite my efforts, I am unable to pinpoint the exact cause of the problem.

Could you please help me identify the issue or suggest how to improve my implementation? Any advice would be greatly appreciated.

Thank you in advance for your help!

Best regards,

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [March 26, 2025, 8:29am UTC](https://getfem.discourse.group/t/help-with-solving-transient-equation-using-getfem-gwfl-approach/70/2 "2025-03-26T08:29:37Z")

</div>

For backward Euler time integration, it should work the way you did it.

Just check the changes found in this version

```python
import numpy as np
import getfem as gf
import os

gf.util_trace_level(1)
gf.util_warning_level(1)

# Paramtres du probleme
L = 1.0 # size of the cube
NX = 40 # mesh szie
NZ = 10 # mesh szie
dt = 2. # time step
Tmax = 40.0 # Temps final

# FitzHugh-Nagumo parametera
Du = 0.01 # Coefficient de diffusion pour u
epsilon = 0.1 # Paramètre epsilon
gamma1 = 0.5 # Paramètre gamma
beta1 = 0.1 # Paramètre beta
a = 0.1 # Seuil d'activation

# mesh
m = gf.Mesh('cartesian', np.linspace(0, L, NX),
                         np.linspace(0, L, NX),
                         np.linspace(0, L/4, NZ))

# FE spaces
mf = gf.MeshFem(m, 1)
mf.set_classical_fem(1) # Élément P1

# integration
mim = gf.MeshIm(m, 4)    

# Model
md = gf.Model("real")

# Variables
md.add_fem_variable("u", mf)
md.add_fem_variable("v", mf)
md.add_fem_data("Previous_u", mf)
md.add_fem_data("Previous_v", mf)
# Initialisation 
# Impulsion initiale au centre
md.set_variable("Previous_u", md.interpolation("-sqr(X(1)-0.5)*sqr(X(2)-0.5)*sqr(X(3)-0.125)", mf))
md.set_variable("Previous_v", md.interpolation("0", mf))

#diffusion terme
md.add_initialized_data("Du", Du)
md.add_initialized_data("dt", dt)

# nonlinear part
md.add_initialized_data("a", a)
md.add_initialized_data("gamma1", gamma1)
md.add_initialized_data("beta1", beta1)
md.add_initialized_data("epsilon", epsilon)

md.add_nonlinear_term(mim, "(u-Previous_u)*Test_u + dt*Du*Grad(u).Grad(Test_u) - dt*(u*(1-u)*(u-a)-v)*Test_u")
md.add_nonlinear_term(mim, "(v-Previous_v)*Test_v - dt*epsilon*(u-gamma1*v+beta1)*Test_v")

if not os.path.exists('Resultf'):
    os.makedirs('Resultf')

# time loop
t = 0
n = 0
while t < Tmax:
    print(f"Temps {t:.2f}")
    md.solve("noisy")
    mf.export_to_vtu(f"Resultf/step_{n:04d}.vtu",
                     mf, md.variable("u"), "u", mf, md.variable("v"), "v")
    md.set_variable("Previous_u", md.variable("u"))
    md.set_variable("Previous_v", md.variable("v"))
    t += dt
    n += 1

print("Simulation terminée. Les résultats sont exportés en VTU.")

```

---

<div class="post-metadata">

**Author:** ![Djoumessi](https://avatars.discourse-cdn.com/v4/letter/d/da6949/32.png) [@Djoumessi](https://getfem.discourse.group/u/Djoumessi)\
**Post date:** [March 27, 2025, 3:50am UTC](https://getfem.discourse.group/t/help-with-solving-transient-equation-using-getfem-gwfl-approach/70/3 "2025-03-27T03:50:28Z")

</div>

Thank very much for your help.  
now it is working.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [March 28, 2025, 8:34am UTC](https://getfem.discourse.group/t/help-with-solving-transient-equation-using-getfem-gwfl-approach/70/4 "2025-03-28T08:34:02Z")

</div>

Nice, I think the main issue was lacking a “minus” after the integration by parts.

Apart from that, please consider the following recommendations:

- Use double quotes for strings, single quote is used for matrix transposition in GWFL.
- Prefer using `model.interpolation()`, which is implemented in GWFL, than the older feature `mesh_fem.eval()` which is implemented with python/numpy.
- Avoid using the old bricks, with GWFL you don’t need bricks.
- Reduce the warnings level with `gf.util_trace_level(1)` and `gf.util_warning_level(1)`, the default is too verbose.
- Use “noisy” in `model.solve()` the default is too silent.
- Avoid using “black\_box” methods like `model.shift_variables_for_time_integration` and `model.add_Newmark_scheme` and similar. All this can be easily implemented in a much more transparent and tidy way with GWFL.
- Prefer the more modern vtu file format than the legacy vtk format for Paraview.
- Learn about partial mesh\_fem and the reduction and extension matrices, it is a very useful concept.

---

<div class="post-metadata">

**Author:** ![Djoumessi](https://avatars.discourse-cdn.com/v4/letter/d/da6949/32.png) [@Djoumessi](https://getfem.discourse.group/u/Djoumessi)\
**Post date:** [March 30, 2025, 7:53pm UTC](https://getfem.discourse.group/t/help-with-solving-transient-equation-using-getfem-gwfl-approach/70/5 "2025-03-30T19:53:39Z")

</div>

Thank very much for your advises. I was able to solve the problem and use that to solve my original problem which was moore complexe than the minimal exemple I’ve provided.  
actually I have some questions if you can helps I would real appreciate.  
I was trying to interpolate a function, like:  
d\*(z\<c)\*(y\<a)\*(x\<b) where a, b, c, d are constants. I’ve tried with `model.interpolation()` and it did not work and my only option was :

```auto
U = mf.eval('1.0*(z < 5)')
md.set_variable('Previous_ui', U)

```

is there any other possibility to write this.

Apart from that, the code is verry slow, is there anywhere to speed-up the code.

saving format vtu or vtk provide too many vtu files, I would like to know if GetFEM has a possibility for hdf5 file for xdmf format.

Thank you very much.

Best.

 ![uLL](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/25b3c8f7d5671aadc50781b68f60317561d663e7.png)

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [March 31, 2025, 5:47am UTC](https://getfem.discourse.group/t/help-with-solving-transient-equation-using-getfem-gwfl-approach/70/6 "2025-03-31T05:47:13Z")

</div>

The syntax in GWFL is `Heaviside(5-X(3))` instead of `z<5`