Help with solving transient equation using GetFEM (GWFL approach)

Djoumessi · March 24, 2025, 8:24pm

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.

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:

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,

Konstantinos.Poulios · March 26, 2025, 8:29am

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

Just check the changes found in this version

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.")

Djoumessi · March 27, 2025, 3:50am

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

Konstantinos.Poulios · March 28, 2025, 8:34am

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.

Djoumessi · March 30, 2025, 7:53pm

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 :

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.

Konstantinos.Poulios · March 31, 2025, 5:47am

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

Topic		Replies	Views
Coordinate naming in MeshFem Questions and issues	5	51	December 11, 2025
Resources for Learning GetFEM++ with Python and GWFL General	12	127	April 1, 2025
Fiber direction import from VTU files for anisotropic strain energy models in GetFEM General	9	54	April 22, 2025
Problem solving transient brick with MUMPS Questions and issues	2	65	January 21, 2025
Implementation of the Gonzalez time scheme for dynamic hyperelasticity Questions and issues	16	126	February 13, 2025

Help with solving transient equation using GetFEM (GWFL approach)

Related topics