Mixed-dimensional branch - Problem with time integration

FEM_TT · February 15, 2023, 8:33am

Hi fenics community,

I use the mixed dimensional branch of @cdaversin and try to solve a time-dependent problem with the following reduced code:

from dolfin import *
from ufl.algorithms import expand_derivatives


# Create mesh
L = 0.20
H = 0.01

mesh = RectangleMesh(Point(0., 0.), Point(L, H), 160, 8, "right")


# Define boundaries and boundary conditions
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
marker.set_all(99)

class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0 + DOLFIN_EPS)


class BoundaryRight(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], L - DOLFIN_EPS)  # H-h


boundary_left = BoundaryLeft()
boundary_left.mark(marker, 1)
boundary_mesh_left = MeshView.create(marker, 1)

boundary_right = BoundaryRight()
boundary_right.mark(marker, 2)

dV = Measure("dx", domain=mesh)
dS_l = Measure("dx", domain=boundary_mesh_left)
dS_r = Measure("ds", domain=mesh, subdomain_data=marker)

tau = Expression(('0.', 't <= tc ? 2*t/tc : 0'), t=0, tc=0.01, degree=1)
nue = Expression(('0', '0'), degree=1)


# Define function space
V_v = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
V_S = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2), symmetry=True)
MX_elem = MixedElement([V_v, V_S])
MX = FunctionSpace(mesh, MX_elem)

LM_left = VectorFunctionSpace(boundary_mesh_left, 'CG', 1, dim=2)

W = MixedFunctionSpace(MX, LM_left)


# Define trial and test functions
w_ = Function(W)
mx, l_left = w_.split()
v, S = split(mx)

dmx, dl_left = TestFunctions(W)
dv, dS = split(dmx)


# Variables of current time step
w_n = Function(W)
mx_n, l_left_n = w_n.split()
v_n, S_n = split(mx_n)


# Time-stepping parameters
T = 0.1
n = 1000
dt = float(T / n)


# Time-stepping
for i in range(n):

    print("Time: ", dt * (i + 1))

    # Update Neumann boundary conditions
    tau.t = float(dt * (i + 1 / 2))

    # Residual
    v_mid = (v + v_n) / 2
    S_mid = (S + S_n) / 2
    l_left_mid = (l_left + l_left_n) / 2

    R0 = dot(dv, (v - v_n) / dt) * dV + \
         inner(grad(dv).T, S_mid) * dV - \
         dot(dv, (l_left - l_left_n) / dt) * dS_l - dot(dv, tau) * dS_r(2)

    R1 = dot(dl_left, (v - v_n) / dt) * dS_l - dot(dl_left, nue) * dS_l

    R2 = inner(dS, (S - S_n) / dt) * dV - inner(dS, grad(v_mid)) * dV

    R = R0 + R1 + R2

    # Solve for next step
    F_blocks = extract_blocks(R)
    Jacobi = []
    for Fi in F_blocks:
        for ui in w_.split():
            d = derivative(Fi, ui)
            d = expand_derivatives(d)
            Jacobi.append(d)

    problem = MixedNonlinearVariationalProblem(F_blocks, w_.split(), [], J=Jacobi)
    solver = MixedNonlinearVariationalSolver(problem)
    solver.solve()

    # Update field
    mx_n.assign(w_.sub(0))
    l_left_n.assign(w_.sub(1))

Unfortunately, the program crashes after a few time iterations and generates a RuntimeError in solver.solve() with the error message:

*** -------------------------------------------------------------------------
*** Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider calling “Function::set_allow_extrapolation(true)” on this Function to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: 68d9ed51bee6723163e7c4fbfffe31e6e3079bca
*** -------------------------------------------------------------------------

It is very strange, because the error does not appear every time and sometimes the code works fine. I wonder why the error occurs randomly.

I appreciate any kind of help to solve my problem

dokken · February 15, 2023, 10:28am

This was most likely fixed in: fenics-project / DOLFIN / commit / 444128a456b9 — Bitbucket

FEM_TT · February 15, 2023, 1:06pm

Hi dokken,

thank you very much for your respond. Is there a docker container including this last commit?
I pulled the docker container mixed dimensional branch of @cdaversin 2 weeks ago.

dokken · February 15, 2023, 1:23pm

See for instance: Package fenics · GitHub

FEM_TT · February 15, 2023, 3:35pm

Hi @dokken,

thank you once again, it seems to work but my code is running much slower now.

dokken · February 15, 2023, 4:55pm

Are you running on a Mac? Then you should use the arm image. If not use the non-arm image

FEM_TT · February 15, 2023, 5:17pm

I am running on windows via pycharm and I use the arm image. I switched to the non-arm image. It is way better now, thank you .

Nonetheless, my original code does not converge now. Did I mess up something with the variable assignment (like in the last lines of my code)? Do I have to consider something special compared to the branch of @cdaversin?

dokken · February 15, 2023, 5:53pm

I dont think so. Could you post an MWE?

FEM_TT · February 16, 2023, 7:17am

Hi @dokken,

for example the following code does not converge at the second time step, when I use the new image.

from dolfin import *
import numpy as np
from ufl.algorithms import expand_derivatives


# Material properties
model = 'plane_stress'
E = 4.5e6
nu = 0.4
rho = 28
Lambda = E * nu / (1 + nu) / (1 - 2 * nu)
G = E / (2.0 * (1.0 + nu))
mu = G
if model == 'plane_stress':
    Lambda = 2 * mu * Lambda / (Lambda + 2 * mu)


# Create mesh
L = 0.20
H = 0.01

mesh = RectangleMesh(Point(0., 0.), Point(L, H), 160, 8, "right")


# Define boundaries and boundary conditions
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
marker.set_all(99)

class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0 + DOLFIN_EPS)


class BoundaryRight(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], L - DOLFIN_EPS) 


boundary_left = BoundaryLeft()
boundary_left.mark(marker, 1)
boundary_mesh_left = MeshView.create(marker, 1)

boundary_right = BoundaryRight()
boundary_right.mark(marker, 2)

dV = Measure("dx", domain=mesh)
dS_l = Measure("dx", domain=boundary_mesh_left)
dS_r = Measure("ds", domain=mesh, subdomain_data=marker)

tau = Expression(('0.', 't <= tc ? 2e3*t/tc : 0'), t=0, tc=0.01, degree=1)
nue = Expression(('0', '0'), degree=1)

# Define function space
V_v = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
V_S = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2), symmetry=True)
V_F = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2))
V_u = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
MX_elem = MixedElement([V_v, V_S, V_F, V_u])
MX = FunctionSpace(mesh, MX_elem)

LM_left = VectorFunctionSpace(boundary_mesh_left, 'CG', 1, dim=2)

W = MixedFunctionSpace(MX, LM_left)


# Define trial and test functions
w_ = Function(W)
mx, l_left = w_.split()
v, S, F, u = split(mx)

dmx, dl_left = TestFunctions(W)
dv, dS, dF, du = split(dmx)


# Variables of current time step
w_n = Function(W)
mx_n, l_left_n = w_n.split()
v_n, S_n, F_n, u_n = split(mx_n)


# Assign initial conditions
MX_F = FunctionSpace(mesh, V_F)
assign((w_.sub(0)).sub(2), project(Identity(MX.mesh().geometry().dim()), MX_F))
assign((w_n.sub(0)).sub(2), project(Identity(MX.mesh().geometry().dim()), MX_F))


# Time-stepping parameters
T = 0.1
n = 1000
dt = float(T / n)


# Time-stepping
for i in range(n):

    print("Time: ", dt * (i + 1))

    # Update Neumann boundary conditions
    tau.t = float(dt * (i + 1 / 2))

    # Residual
    v_mid = (v + v_n) / 2
    S_mid = (S + S_n) / 2
    F_mid = (F + F_n) / 2
    l_left_mid = (l_left + l_left_n) / 2

    R0 = dot(dv, rho * (v - v_n) / dt) * dV + \
         inner(dot(grad(dv).T, F_mid), S_mid) * dV - \
         dot(dv, (l_left - l_left_n) / dt) * dS_l - dot(dv, tau) * dS_r(2)

    R1 = dot(dl_left, (v - v_n) / dt) * dS_l - dot(dl_left, nue) * dS_l

    R2 = inner(dS, -(Lambda / (2 * mu * (3 * Lambda + 2 * mu))) * tr((S - S_n) / dt) * Identity(2) + 1 / (2 * mu) * (S - S_n) / dt) * dV - inner(dS, dot(F_mid.T, grad(v_mid))) * dV


    R3 = inner(dF, (F - F_n) / dt) * dV - \
     inner(dF, grad(v_mid)) * dV

    R4 = dot(du, (u - u_n) / dt) * dV - \
     dot(du, v_mid) * dV

    R = R0 + R1 + R2 + R3 + R4

    # Solve for next step
    F_blocks = extract_blocks(R)
    Jacobi = []
    for Fi in F_blocks:
        for ui in w_.split():
            d = derivative(Fi, ui)
            d = expand_derivatives(d)
            Jacobi.append(d)

    problem = MixedNonlinearVariationalProblem(F_blocks, w_.split(), [], J=Jacobi)
    solver = MixedNonlinearVariationalSolver(problem)
    solver.solve()

    # Update field
    mx_n.assign(w_.sub(0))
    l_left_n.assign(w_.sub(1))

It is strange and I have probally mess up something with the variable assignment or the split function?

dokken · February 16, 2023, 7:46am

This is far from and MWE (minimal working example), and I do not have time to parse all this code.

One thing that stands out is:
Why do you redefine the whole variational form inside the loop?

FEM_TT:

    v_mid = (v + v_n) / 2
    S_mid = (S + S_n) / 2
    F_mid = (F + F_n) / 2
    l_left_mid = (l_left + l_left_n) / 2

    R0 = dot(dv, rho * (v - v_n) / dt) * dV + \
         inner(dot(grad(dv).T, F_mid), S_mid) * dV - \
         dot(dv, (l_left - l_left_n) / dt) * dS_l - dot(dv, tau) * dS_r(2)

    R1 = dot(dl_left, (v - v_n) / dt) * dS_l - dot(dl_left, nue) * dS_l

    R2 = inner(dS, -(Lambda / (2 * mu * (3 * Lambda + 2 * mu))) * tr((S - S_n) / dt) * Identity(2) + 1 / (2 * mu) * (S - S_n) / dt) * dV - inner(dS, dot(F_mid.T, grad(v_mid))) * dV


    R3 = inner(dF, (F - F_n) / dt) * dV - \
     inner(dF, grad(v_mid)) * dV

    R4 = dot(du, (u - u_n) / dt) * dV - \
     dot(du, v_mid) * dV

    R = R0 + R1 + R2 + R3 + R4

    # Solve for next step
    F_blocks = extract_blocks(R)
    Jacobi = []
    for Fi in F_blocks:
        for ui in w_.split():
            d = derivative(Fi, ui)
            d = expand_derivatives(d)
            Jacobi.append(d)

    problem = MixedNonlinearVariationalProblem(F_blocks, w_.split(), [], J=Jacobi)
    solver = MixedNonlinearVariationalSolver(problem)

and Im not certain about your usage of split,

I believe that should be split(w). This holds in all the different places you use split

FEM_TT · February 16, 2023, 8:30am

Hi @dokken,

I am sorry, I tried to create a MWE.

Using the following code

from dolfin import *

# Create mesh
mesh = RectangleMesh(Point(0., 0.), Point(0.2, 0.01), 160, 8, "right")


# Define boundaries
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
marker.set_all(99)

class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0 + DOLFIN_EPS)

boundary_left = BoundaryLeft()
boundary_left.mark(marker, 1)
boundary_mesh_left = MeshView.create(marker, 1)

dV = Measure("dx", domain=mesh)
dS_l = Measure("dx", domain=boundary_mesh_left)


# Define function space
V_v = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
V_S = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2), symmetry=True)
MX_elem = MixedElement([V_v, V_S])
MX = FunctionSpace(mesh, MX_elem)

LM_left = VectorFunctionSpace(boundary_mesh_left, 'CG', 1, dim=2)

W = MixedFunctionSpace(MX, LM_left)


# Define trial and test functions
w_ = Function(W)
mx, l_left = w_.split()
v, S = split(mx)

dmx, dl_left = TestFunctions(W)
dv, dS = split(dmx)


# Weak form
tau = Expression(('0.', '20'), degree=1)

R0 = dot(dv, v) * dV + inner(grad(dv).T , S) * dV - \
     dot(dv, l_left) * dS_l - dot(dv, tau) * dV

R1 = dot(dl_left, v) * dS_l

R2 = inner(dS, tr(S) * Identity(2) + S) * dV - \
     inner(dS, grad(v)) * dV

R = R0 + R1 + R2


# Solution
solve(R==0, w_)

works well using the branch of @cdaversin. But the solver does not converge in case of the new image.

IngeborgGjerde · February 21, 2023, 12:23pm

Hi,

I tested the MWE posted above and I also cannot get it to work on ghcr.io/scientificcomputing/fenics:2023-02-17 or on ceciledc/fenics_mixed_dimensional:16-01-23. As FEM_TT says the Newton solver fails to converge. I will post a bug report on bitucket.

It does work and give the correct result using e.g. ceciledc/fenics_mixed_dimensional:21-06-22 - but then your original problem remains wrt the program crashing

You could assemble the code using fenics_ii instead:

from xii import *
from xii.nonlin.jacobian import block_jacobian
from fenics import *

# Create mesh
mesh = RectangleMesh(Point(0., 0.), Point(0.2, 0.01), 160, 8, "right")

# Define boundaries
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
marker.set_all(99)

class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0 + DOLFIN_EPS)

boundary_left = BoundaryLeft()
boundary_left.mark(marker, 1)
boundary_mesh_left = EmbeddedMesh(marker, 1)

dV = Measure("dx", domain=mesh)
dS_l = Measure("dx", domain=boundary_mesh_left)


# Define function space
V_v = VectorFunctionSpace(mesh, 'CG', 1, 2)
V_S = TensorElement('DG', mesh.ufl_cell(), 0, shape=(2, 2), symmetry=True)
V_S = FunctionSpace(mesh, V_S)

LM_left = VectorFunctionSpace(boundary_mesh_left, 'CG', 1, dim=2)

W = [V_v, V_S, LM_left] # list of function spaces

# Define trial and test functions
w = ii_Function(W) 
v, S, l_left = w # w is a list of functions

dv, dS, dl_left = map(TestFunction, W)

Tdv = Trace(dv, boundary_mesh_left)
Tv = Trace(v, boundary_mesh_left)

# Weak form
tau = Expression(('0.', '20'), degree=1)
tau = interpolate(tau, V_v) #fenics_ii later needs to know the "type" of tau, that is why we interpolate it

R0 = dot(dv, v) * dV + inner(grad(dv).T , S) * dV + \
     - dot(Tdv, l_left) * dS_l \
     - dot(dv, tau) * dV

R1 = dot(dl_left, Tv) * dS_l

R2 = inner(dS, tr(S) * Identity(2) + S) * dV - \
     inner(dS, grad(v)) * dV

F = [R0, R2, R1]

# Solution
dF = block_jacobian(F, w)

# Newton
eps = 1.0
tol = 1.0E-10
niter = 0
maxiter = 125

dw = ii_Function(W)
while eps > tol and niter < maxiter:
    niter += 1

    A, b = list(map(ii_assemble, (dF, F)))
    A, b = list(map(ii_convert, (A, b)))

    solve(A, dw.vector(), b)

    eps = sqrt(sum(x.norm('l2')**2 for x in dw.vectors()))

    print('\t%d |du| = %.6E |A|= %.6E |b| = %.6E' % (niter, eps, A.norm('linf'), b.norm('l2')))

    for i in range(len(W)):
        w[i].vector().axpy(-1, dw[i].vector())

You can install fenics_ii using pip (as described in the readme, feel free to message me if you need assistance) or run on the docker image mirok/emi-book-fem

FEM_TT · February 23, 2023, 10:14am

Hi @IngeborgGjerde,

thank you very much, that is very kind.
Unfortunately, I can not find the modul “xii” using the docker image mirok/emi-book-fem:latest. Do you know, what might be the problem?

IngeborgGjerde · February 26, 2023, 2:35pm

Aha, the mirok/emi-book-fem image is a bit old and uses python2.

Can you try this image: ghcr.io/ingeborggjerde/graphnics:v0.8.0 ?

Sorry for the late reply, I am travelling at the moment. I’ll update my post in a few days with pip installation instructions

FEM_TT · February 27, 2023, 10:37am

Hi @IngeborgGjerde,

thank you again, your image seems to work perfectly. Now, I will try to solve my original problem and I will let you know, if it works . I wish you a nice trip.

Topic		Replies	Views
Mixed-dimensional branch - Solving optimal control problem	2	899	December 16, 2019
[Mixed-dimensional Branch] Coupling Subdomain function with whole-domain function Errors	1	601	May 3, 2022
Different result when using fenics_mixed_dimensional branch Errors	1	345	July 7, 2021
Mixed-dimensional branch - Error: Unable to set given (local) rows to identity matrix Legacy DOLFIN	2	124	May 18, 2024
timeNS L2,H1 error	2	579	April 21, 2022

Mixed-dimensional branch - Problem with time integration

Related topics