Explicit time solver

General
1

Hi all,

This is may be a silly question. I wanted to compare the performance of an implicit solver I built in dolfinx for the linear elastodynamic equations ( \rho \frac{\partial^2 u}{\partial t^2} = \nabla \cdot \sigma ), to an explicit solver e.g using velocity verlet to update time. This completely removes the need to solve any linear system or non linear system as the new acceleration is computed from \rho a_{n+1}=\nabla\cdot\sigma(u_{n+1}) and then displacement is updated through a linear combination of the old quantities.

I haven’t seen much information on explicit time solvers in the FEniCS project, so I was wondering whether I can do it in dolfinx and still utilise nice features like MPI parallelisation? Or is this platform not really built to solve problems in this way?

Any information/ help would be greatly appreciated.

Thanks :slight_smile:

2

Hi, you might have a look on the following code:

from dolfinx.fem.petsc import assemble_vector
from dolfinx.fem import Function, set_bc, form
from ufl import action

def petsc_div(numerateur, denominateur, output):
    """ 
    Pointwise division between two vectors using PETSc. Le vecteur
    output est rempli avec le résultat x/y

    numerateur : PETScVector
    denominateur : PETScVector
    output : PETScVector
    """
    output.pointwiseDivide(numerateur, denominateur)
    
def dt_update(x, dot_x, dt, new_vec = False):
    """
    Mise jour explicite de x en utilisant sa dérivée, Euler-explicite

    Parameters
    ----------
    x : Function, fonction Ă  mettre Ă  jour.
    dot_x : Function, dérivée temporelle de x.
    dt : Float, pas de temps temporel.
    """
    if new_vec:
        u = x.copy()
        u.vector.axpy(dt, dot_x.vector)
        return u
    else:
        x.vector.axpy(dt, dot_x.vector)
        return

class ExplicitDisplacementSolver:
    def __init__(self, u, v, a, dt, m_form, residual, bcs):
        """
        Initilise le solveur explicit de calcul du déplacement

        Parameters
        ----------
        u : Function, champ de déplacement
        v : Function, champ de vitesse
        a : Function, champ d'accélération
        dt : float, pas de temps
        m_form : Form, forme biilinéaire de masse
        residual : Form, résidu du problème mécanique: a(u,v) - L(v)
        bcs : DirichletBC, conditions aux limites en déplacement
        """
        self.u = u
        self.v = v
        self.a = a
        self.dt = dt
        self.bcs = bcs
        self.set_explicit_function(residual, m_form)
        
    def set_explicit_function(self, residual_form, m_form):
        """
        DĂ©finition du vecteur de masse, issu de la condensation 
        de la matrice de masse
        
        Parameters
        ----------
        residual_form: Form, résidu = forme linéaire a(u,v)-l(v)
        m_form : Form, forme biilinéaire de masse.
        """
        u1 = Function(self.u.function_space)
        u1.vector.array[:] = 1.
        self.diag_M = assemble_vector(form(action(m_form, u1)))
        self.local_res = form(- residual_form)
        self.compute_acceleration_speed()

    def compute_acceleration_speed(self):
        """
        Calcul le résidu puis l'utilise afin de 
        calculer la vitesse et l'accélération
        """
        res = assemble_vector(self.local_res)
        petsc_div(res, self.diag_M, self.a.vector)
        set_bc(self.a.vector, self.bcs.bcs_homog)
        dt_update(self.v, self.a, self.dt)
        res.destroy()
            
    def u_solve(self):
        """
        Calcl du déplacement sur une iteration explicite:
        -actualisation des déplacements
        -imposition des CLs en déplacements
        -calcul de l'acceleration et de la vitesse
        """      
        dt_update(self.u, self.v, self.dt)
        set_bc(self.u.vector, self.bcs.bcs)
        self.compute_acceleration_speed()

I hope it can help,
Paul

1 Like
3

Thanks Paul, this was really useful. I have a working solver for my setup. Is it simple to parallelise this solver - simply running on more than one core seems to go wrong for me?

1 Like
4

I should expand…

I assumed the issue was I needed to update the ghost data after the petsc_div with self.a.x.scatter_forward() in the compute_acceleration_speed function. But it doesn’t seem to work

Any help would be appreciated :slight_smile:

5

@BouteillerP, the following MWE gives different results for parallel and series computation:

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, default_scalar_type, mesh
from dolfinx.io import XDMFFile
import ufl
import dolfinx.fem.petsc

import time
start_time = time.time()

#============= Create mesh =============
L = 1.0  # Length of the rod
nx, ny, nz = 10, 10, 100  # Number of elements in x, y, and z directions
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array(
    [0.1, 0.1, L])], [nx, ny, nz], cell_type=mesh.CellType.hexahedron)

#============= Create function spaces =============
V = fem.VectorFunctionSpace(domain, ("CG", 1))

fdim = domain.topology.dim - 1

#============= Physical constants =============
E = 200e9
nu = 0.27
mu = E / (2.0 * (1 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
rho = 7990
Amp=1e-8    # Amplitude of excitation

#============= Time integration parameters =============
T = 5.0e-4
NTsteps = 2000
save_every = 20
dt = fem.Constant(domain, T / NTsteps)

#============= Time dep. boundary conditions =============
def point1(x):
    return np.logical_and(np.logical_and(np.isclose(x[2], 0.0), np.isclose(x[1], 0.05)), np.isclose(x[0], 0.0))
def point2(x):
    return np.logical_and(np.logical_and(np.isclose(x[2], 0.0), np.isclose(x[1], 0.05)), np.isclose(x[0], 0.1))

def pert(t):
    return Amp*np.sin(2*np.pi*t/T)

def input1(t):
    return np.array([0.0,pert(t),0.0], dtype=default_scalar_type)
def input2(t):
    return np.array([0.0,-pert(t),0.0], dtype=default_scalar_type)

bc_fun1=fem.Constant(domain, PETSc.ScalarType((0,0,0)))
bc_fun2=fem.Constant(domain, PETSc.ScalarType((0,0,0)))

point1_dofs = fem.locate_dofs_geometrical(V, point1) 
point2_dofs = fem.locate_dofs_geometrical(V, point2) 
 
bc1 = fem.dirichletbc(bc_fun1, point1_dofs, V)
bc2 = fem.dirichletbc(bc_fun2, point2_dofs, V)

allbcs=[bc1,bc2]

#============= Functions and fields =============
u_ = ufl.TestFunction(V)
u = fem.Function(V, name='Displacement')
v = fem.Function(V)
a = fem.Function(V)
a_ = ufl.TrialFunction(V)

#============= Measures =============
dx = ufl.Measure("dx", domain=domain)

#============= Stress and bilinear forms =============
def sigma(r):
    return lmbda*ufl.nabla_div(r)*ufl.Identity(len(r)) + 2*mu*ufl.sym(ufl.grad(r))

def m(a, u_):
    return rho*ufl.inner(a, u_)*dx

def k(u, u_):
    return ufl.inner(sigma(u), ufl.grad(u_))*dx

#============= Explicit solver functions setup =============
def petsc_div(numerateur, denominateur, output):
    """ 
    Pointwise division between two vectors using PETSc.

    numerateur : PETScVector
    denominateur : PETScVector
    output : PETScVector
    """
    output.pointwiseDivide(numerateur, denominateur)

#============= Residual to solve =============
residual =  k(u, u_)
m_form = m(a_,u_)

#============= Create file =============
xdmf_file=XDMFFile(domain.comm, "output_xdmfs/testnocores.xdmf", "w")

xdmf_file.write_mesh(domain)
xdmf_file.write_function(u, 0)

#============= Time loop =============
tt=fem.Constant(domain, 0.)

u1 = fem.Function(V)
u1.vector.array[:] = 1.
diag_M = fem.petsc.assemble_vector(fem.form(ufl.action(m_form, u1)))
diag_M.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
diag_M.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
local_res = fem.form(-residual)

for i in range(NTsteps):
    tt.value += dt.value
    
    bc_fun1.value = input1(tt.value)
    bc_fun2.value = input2(tt.value)
    
    #------------ Explicit Time stepping --------------
    u.x.array[:]=u.x.array+ float(dt)*v.x.array
    fem.set_bc(u.vector, allbcs)
    u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    local_res = fem.form(-residual)
    res = fem.petsc.assemble_vector(local_res)
    res.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    res.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    petsc_div(res, diag_M, a.vector)
    a.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    v.x.array[:]=v.x.array+ float(dt)*a.x.array
    v.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    res.destroy()

    if MPI.COMM_WORLD.rank==0:
        print(f"Time: {tt.value:.5e}")

    #------------- Write variables to XDMF file -------------
    if (i+1) % save_every == 0:
        xdmf_file.write_function(u, tt.value)

xdmf_file.close()

if MPI.COMM_WORLD.rank==0:
    print("--- %s seconds runtime ---" % (time.time() - start_time))
6

I would strongly suggest to reduce this to a minimal example.
For instance:

#!/usr/bin/env python3

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, mesh
from dolfinx.io import XDMFFile
import dolfinx.fem.petsc


L = 1.0  # Length of the rod
nx, ny, nz = 3, 3, 3  # Number of elements in x, y, and z directions
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array(
    [0.1, 0.1, L])], [nx, ny, nz], cell_type=mesh.CellType.hexahedron)


def petsc_div(numerateur, denominateur, output):
    """ 
    Pointwise division between two vectors using PETSc.

    numerateur : PETScVector
    denominateur : PETScVector
    output : PETScVector
    """
    output.pointwiseDivide(numerateur, denominateur)


xdmf_file = XDMFFile(MPI.COMM_WORLD, "test_mpi.xdmf", "w")
xdmf_file.write_mesh(domain)

V = fem.VectorFunctionSpace(domain, ("CG", 1))
u = dolfinx.fem.Function(V)
a = dolfinx.fem.Function(V)
res = dolfinx.fem.Function(V)
for i in range(1, 10):
    a.interpolate(lambda x: (x[0], x[1], x[2]))
    u.interpolate(lambda x: (np.sin(x[0]), np.sin(x[1]), np.sin(x[2])))
    # ------------ Explicit Time stepping --------------
    petsc_div(a.vector, res.vector, u.vector)
    u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                         mode=PETSc.ScatterMode.FORWARD)
    u.x.array[:] = u.x.array + float(i)*a.x.array
    u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                         mode=PETSc.ScatterMode.FORWARD)

    xdmf_file.write_function(u, i)

xdmf_file.close()

does not reproduce the issue in serial and parallel, so I would suggest that you reduce the problem to the smallest amount of code that reproduces the parallel problem.

7

Apologies Dokken, I will bear this in mind for next time. However, I have found the source of issue…

When setting up the linear mass form (pre PETsc division),

u1 = fem.Function(V)
u1.vector.array[:] = 1.
diag_M = fem.petsc.assemble_vector(fem.form(ufl.action(m_form, u1)))

should be replaced with

u1 = fem.Function(V)
diag_m_form = fem.form(ufl.action(m_form, u1))
u1.x.array[:] = 1.
diag_M = fem.petsc.assemble_vector(diag_m_form)

I’m not entirely sure why “.x” versus “.vector” is necessary for consistency in parallel.

8

dolfinx.fem.Function.x is the DOLFINx implementation of a partitioned array. The array you access through dolfinx.fem.Function.x.array includes all owned entries + ghosted entries.

dolfinx.fem.Function.vector returns a PETSc representation of the same array. Their access to non-local entries are more complicated, and you would have to explicitly set ghost entries, or use scatter forward after setting data to the array.

2 Likes