Initial conditions for mixed elements with vector and scalar elements

General
1

Hello Fenics community,

I am trying to implement an RD system where I also have to update the mesh coordinates. When I try to assign the initial conditions for the mesh coordinate function.

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

def x_expr(x):
    """Expression for the exact velocity solution to Kovasznay flow"""
    return np.vstack((x[0],x[1]))

x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())

normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)

#X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)
X0.sub(0).interpolate(x_expr)

X0.x.scatter_forward()

I obtain strange values if I plot as a vector (It should look like outward vectors).

image
image838×487 47.2 KB

I have already implemented it in different ways, but I cannot find the error.

Here is my MWE

# +
import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math
import matplotlib.pyplot as plt

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot

# Save all logging to file
log.set_output_file("log.txt")
# - 

# Next, various model parameters are defined:

dt = 5.0e-02  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.3)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)

    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.
    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))

n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

signx = ufl.sign(ufl.dot(n, r))
xc = signx*n

V1 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)
V2 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)

ME2 = functionspace(mesh, mixed_element([V1, V2],gdim=2))

# Trial and test functions of the space `ME` are now defined:

q, v = ufl.TestFunctions(ME)

Xi, phi = ufl.TestFunctions(ME2)

# ```{index} split functions
# ```

# +
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

# X = Function(ME2)  # current solution
X0 = Function(ME2)  # solution from previous converged step
X0.name = "X0"

(Xh, H) = ufl.TrialFunctions(ME2)  # current solution
#(Xh0, H0) = ufl.TrialFunctions(ME2)  # solution from previous converged step

# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)

# # Split mixed functions
# Xh, H = ufl.split(X)
Xh0 = X0.sub(0)
H0 = X0.sub(1)

au = 0.9
bu = 1.1
av = 0.8
bv = 0.95

# Interpolate initial condition

u.sub(0).interpolate(lambda x: (bu-au)*np.random.rand(x.shape[1]) +au)
u.sub(1).interpolate(lambda x: (bv-av)*np.random.rand(x.shape[1]) +av)

u.x.scatter_forward()

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

def x_expr(x):
    """Expression for the exact velocity solution to Kovasznay flow"""
    return np.vstack((x[0],x[1]))

x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())

normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)

#X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)
X0.sub(0).interpolate(x_expr)

X0.x.scatter_forward()

print (f"X array  {X0.sub(0).x.array[:]}")

u_vec = ufl.as_vector([u1, u2])
kp = ufl.as_vector([1.5 , 0])

# Weak statement of the equations

k = dt
d = 10
gamma = 29
a = 0.1
b = 0.9
tol = 1e-4

d1 = 1.0
d2 = d

ks = 25
kb = 2

F = ((u1 - u10) / k)*q*dx + d1*inner(grad(u1), grad(q))*dx\
   -(gamma*(u1*u1*u2-u1+a))*q*dx \
  + ((u2 - u20) / k)*v*dx + d2*inner(grad(u2), grad(v))*dx\
  -(gamma*(-u1*u1*u2+b))*v*dx 

ax = dot((Xh / k),n_oriented)*phi*dx  + kb*inner(grad(H),grad(phi))*dx  - 1/2 *kb* H0**2  * H *phi*dx + ks*H*phi*dx\
     + dot(H*n_oriented,Xi)*dx - inner(grad(Xh),grad(Xi))*dx

Lx = dot(kp,u_vec)*phi*dx + dot((Xh0 / k),n_oriented)*phi*dx


# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver

ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

# X_h = Function(ME2)
# X_h.name = "X_h"

problemX = LinearProblem(ax, Lx)

# Step in time
t = 0.0

# Output file
name = "mesh_at_t"+str(t)
mesh.name = name

# Get the sub-space for c and the corresponding dofs 

normal_vec.interpolate(normal_expr)

V0, dofs = ME.sub(0).collapse()
W02, dofs2 = ME2.sub(0).collapse()

u1 = u.sub(0)
u2 = u.sub(1)

Xh0 = X0.sub(0)
H0 = X0.sub(1)

X_coord = X0.x.array[dofs2].reshape(len(mesh.geometry.x[:,1]),2)

file1 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/outputu1.xdmf", "w")
file3 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/normal.xdmf", "w")
file4 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/H.xdmf", "w")

file1.write_mesh(mesh)
file3.write_mesh(mesh)
file4.write_mesh(mesh)

file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file4.write_function(X0.sub(0), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

norm_stop = 0.1e-6

fileL2 = open('L2.dat', 'w') 

np.savetxt("coordinates.dat", mesh.geometry.x, delimiter='\t', header='x\ty', comments='')

z = 0
inc = 0 
impresion = 1

u0.x.array[:] = u.x.array
#print (f"MESH array {mesh.geometry.x[:,]}")

t += dt

X_h = problemX.solve()
r = solver.solve(u)
inc += 1

print(f"Step {int(t/dt)}: num iterations: {r[0]}")

l2_norm = np.linalg.norm(u.x.array[dofs]-u0.x.array[dofs])/dt
#print (f"L2_norm = {l2_norm}")

fileL2.write(f"{t}\t{l2_norm}\n")
u0.x.array[:] = u.x.array
X0.x.array[:] = X_h.x.array

num=inc/impresion-1 


X_coord = X0.x.array[dofs2].reshape(len(mesh.geometry.x[:,1]),2)
DeltaX = mesh.geometry.x[:,[0,1]]-X_coord

#print (f"Delta X array  {DeltaX}")
mesh.geometry.x[:,0] = X_coord[:,0]
mesh.geometry.x[:,1] = X_coord[:,1]

normal_vec.interpolate(normal_expr)

# Save solution to file (VTK)

name = "mesh_at_t"+str(t)
mesh.name = name

file1.write_mesh(mesh)
file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file3.write_mesh(mesh)
file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file4.write_mesh(mesh)
file4.write_function(X0.sub(0), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
z += 1

file1.close()
file3.close()
file4.close()

fileL2.close()

print(l2_norm)
#print(H.x.array)

Does anybody know what could be the error?

Thank you.

2

Try calling X0.sub(0).collapse().

Please note that this is quite far away from a minimal example, as you have provided 4 output files, without specifying which one we should look at, and very complicated for for the mesh reading (Why write it to file if you are reading it in again).

There is also tons of code not relating to the normal vector, but to your full problem.
This makes it extremely hard for anyone to give sensible advice.

3

Sorry about that.

I was trying to focus on the initial conditions, where I think the error is, but I will submit my code again. Actually I am only interested in the function ‘u1’ but I can obtain correct results if I can make that the mesh evolve properly

import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot

# Save all logging to file
log.set_output_file("log.txt")
# -

# Next, various model parameters are defined:

dt = 5.0e-03  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.3)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)

    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.
    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))

n = ufl.CellNormal(mesh)
xc = ufl.SpatialCoordinate(mesh)
r = xc/ufl.sqrt(ufl.dot(xc, xc))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

V1 = element("CG", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)
V2 = element("CG", mesh.basix_cell(), 2, gdim=2)

ME2 = functionspace(mesh, mixed_element([V1, V2],gdim=2))

# Trial and test functions of the space `ME` are now defined:

q, v = ufl.TestFunctions(ME)

Xi, phi = ufl.TestFunctions(ME2)

# ```{index} split functions
# ```

# +
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

X0 = Function(ME2)  # solution from previous converged step

(Xh, H) = ufl.TrialFunctions(ME2)  # current solution

# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)

# # Split mixed functions
Xh0, H0 = ufl.split(X0)

au = 0.9
bu = 1.1
av = 0.8
bv = 0.95

# Interpolate initial condition

u.sub(0).interpolate(lambda x: (bu-au)*np.random.rand(x.shape[1]) +au)
u.sub(1).interpolate(lambda x: (bv-av)*np.random.rand(x.shape[1]) +av)

u.x.scatter_forward()

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

def x_exp(x):
    """Expression for coordinates"""
    return np.vstack((x[0],x[1]))

H_exp = Expression(div(sign*n), V02.element.interpolation_points())

X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)

X0.x.scatter_forward()

u_vec = ufl.as_vector([u1, u2])
kp = ufl.as_vector([1.5 , 0])

# Weak statement of the equations

k = dt
d = 10
gamma = 29
a = 0.1
b = 0.9
tol = 1e-4

d1 = 1.0
d2 = d

ks = 25
kb = 2

F = ((u1 - u10) / k)*q*dx + d1*inner(grad(u1), grad(q))*dx\
   -(gamma*(u1*u1*u2-u1+a))*q*dx \
  + ((u2 - u20) / k)*v*dx + d2*inner(grad(u2), grad(v))*dx\
  -(gamma*(-u1*u1*u2+b))*v*dx 

ax = dot((Xh / k),n_oriented)*phi*dx  + kb*inner(grad(H),grad(phi))*dx  - 1/2 *kb* H0**2  * H *phi*dx + ks*H*phi*dx\
     + dot(H*n_oriented,Xi)*dx - inner(grad(Xh),grad(Xi))*dx

Lx = dot(kp,u_vec)*phi*dx + dot((Xh0 / k),n_oriented)*phi*dx

# Create nonlinear problem and Newton solver

problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver

ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

#Linear problem

problemX = LinearProblem(ax, Lx)

# Step in time
t = 0.0

# Output file
name = "mesh_at_t"+str(t)
mesh.name = name

#  Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
    T = 3 * dt
else:
    T = 10 * dt

# Get the sub-space for u1 and X and the corresponding dofs 


V0, dofs = ME.sub(0).collapse()
W02, dofs2 = ME2.sub(0).collapse()

u1 = u.sub(0)
u2 = u.sub(1)

Xh0 = X0.sub(0)
H0 = X0.sub(1)

file1 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/outputu1.xdmf", "w")
file2 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/H.xdmf", "w")

file1.write_mesh(mesh)
file2.write_mesh(mesh)

file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file2.write_function(X0.sub(1), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

norm_stop = 0.1e-6

fileL2 = open('L2.dat', 'w') 

z = 0
inc = 0 
impresion = 1

u0.x.array[:] = u.x.array

while (t < T):
    t += dt

    X_h = problemX.solve()
    r = solver.solve(u)
    inc += 1

    print(f"Step {int(t/dt)}: num iterations: {r[0]}")

    l2_norm = np.linalg.norm(u.x.array[dofs]-u0.x.array[dofs])/dt

    fileL2.write(f"{t}\t{l2_norm}\n")
    u0.x.array[:] = u.x.array
    X0.x.array[:] = X_h.x.array

    num=inc/impresion-1 

    X_coord = X0.x.array[dofs2].reshape(len(mesh.geometry.x[:,1]),2)

    # Move the mesh

    mesh.geometry.x[:,0] = X_coord[:,0]
    mesh.geometry.x[:,1] = X_coord[:,1]

    if (int(num)==z) or (int(num-z)==0):
        us, vs = ufl.split(u)

        # Save solution to file (VTK)

        name = "mesh_at_t"+str(t)
        mesh.name = name

        file1.write_mesh(mesh)
        file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
        file2.write_mesh(mesh)
        file2.write_function(H0, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

        z += 1
    
file1.close()

fileL2.close()

print(l2_norm)

With the system for moving the mesh, the mesh behaves poorly from the first step. This is why I think I may not be imposing the initial conditions correctly, or perhaps the error could lie in assigning the new mesh coordinates. I think that the weak formulation is good (Xh and Xh0 functions).

image
image923×338 56.2 KB

Thank you,

4

Please the reduce the code to applying the initial condition, removing all other code.

5

Here is the code:

import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot

# Save all logging to file
log.set_output_file("log.txt")
# -

# Next, various model parameters are defined:

dt = 5.0e-03  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.4)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)

    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.
    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))

n = ufl.CellNormal(mesh)
xc = ufl.SpatialCoordinate(mesh)
r = xc/ufl.sqrt(ufl.dot(xc, xc))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

V1 = element("CG", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)
V2 = element("CG", mesh.basix_cell(), 2, gdim=2)

ME2 = functionspace(mesh, mixed_element([V1, V2],gdim=2))

# Trial and test functions of the space `ME` are now defined:

q, v = ufl.TestFunctions(ME)

Xi, phi = ufl.TestFunctions(ME2)

# ```{index} split functions
# ```

# +
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

X0 = Function(ME2)  # solution from previous converged step

(Xh, H) = ufl.TrialFunctions(ME2)  # current solution

# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)

W02, dofs2 = ME2.sub(0).collapse()

# # Split mixed functions
Xh0, H0 = ufl.split(X0)

au = 0.9
bu = 1.1
av = 0.8
bv = 0.95

# Interpolate initial condition

u.sub(0).interpolate(lambda x: (bu-au)*np.random.rand(x.shape[1]) +au)
u.sub(1).interpolate(lambda x: (bv-av)*np.random.rand(x.shape[1]) +av)

u.x.scatter_forward()

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

def x_exp(x):
    """Expression for coordinates"""
    return np.vstack((x[0],x[1]))

H_exp = Expression(div(n_oriented), V02.element.interpolation_points())

X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)

X0.x.scatter_forward()

X_coord = X0.x.array[dofs2].reshape(len(mesh.geometry.x[:,1]),2)
# Move the mesh

mesh.geometry.x[:,0] = X_coord[:,0]
mesh.geometry.x[:,1] = X_coord[:,1]


u1 = u.sub(0)
u2 = u.sub(1)

Xh0 = X0.sub(0)
H0 = X0.sub(1)

# Step in time
t = 0.0

# Output file
name = "mesh_at_t"+str(t)
mesh.name = name

file1 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/outputu1.xdmf", "w")
file2 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/H.xdmf", "w")

file1.write_mesh(mesh)
file2.write_mesh(mesh)

file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file2.write_function(H0, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

I think that the problem could be in the line X0.sub(0).interpolate(x_exp), if you see the output file ‘outputu1.xdmf’ you will see something like this:

image

Thank you in advance

7

@dokken, I’ve just realized what the problem is: I hadn’t taken into account the order of the nodes for second-order elements. When I use first-order elements, my code works well. However, when I use second-order elements, this line of the code performs poorly:

Xh0.x.array[:] = Xn.x.array[:]
H0.x.array[:] = Hn.x.array[:]

num=inc/impresion-1 

X_coord2 = Xh0.x.array.reshape(len(mesh.geometry.x[:,1]),2)
xf2 = Xh0.x.array[:]

mesh.geometry.x[:,0] = X_coord2[:,0]
mesh.geometry.x[:,1] = X_coord2[:,1]

Now, I cannot understand how the nodes are arranged in the mesh and how they are arranged in the elements. I found some information about it, but nothing on how I can map this in my case. Is there a way to relate these different arrays?

Here is the complete code if you want to check

# +
import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot

# Save all logging to file
log.set_output_file("log.txt")
# -

# Next, various model parameters are defined:

dt = 5.0e-05  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.4)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)

    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.
    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))

n = ufl.CellNormal(mesh)
xc = ufl.SpatialCoordinate(mesh)
r = xc/ufl.sqrt(ufl.dot(xc, xc))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)

ME2 = functionspace(mesh, mixed_element([P2, P1],gdim=2))

V, Q = functionspace(mesh, P2), functionspace(mesh, P1)


# Trial and test functions of the space `ME` are now defined:

q, v = ufl.TestFunctions(ME)

Xi, phi = ufl.TestFunctions(ME2)

# ```{index} split functions
# ```

# +
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

Xh0 = Function(V)  # solution from previous converged step
H0 = Function(Q)

(Xh, H) = ufl.TrialFunctions(ME2)  # current solution

# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)

au = 0.9
bu = 1.1
av = 0.8
bv = 0.95

# Interpolate initial condition

u.sub(0).interpolate(lambda x: (bu-au)*np.random.rand(x.shape[1]) +au)
u.sub(1).interpolate(lambda x: (bv-av)*np.random.rand(x.shape[1]) +av)

u.x.scatter_forward()

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

x_exp = Expression(xc,V.element.interpolation_points())
H_exp = Expression(div(sign*n),Q.element.interpolation_points())

Xh0.interpolate(x_exp)
H0.interpolate(H_exp)

Xh0.x.scatter_forward()
H0.x.scatter_forward()

u_vec = ufl.as_vector([u1, u2])
kp = ufl.as_vector([1.5 , 0])

# Weak statement of the equations

k = dt
d = 10
gamma = 29
a = 0.1
b = 0.9
tol = 1e-4

d1 = 1.0
d2 = d

ks = 25
kb = 2

F = ((u1 - u10) / k)*q*dx + d1*inner(grad(u1), grad(q))*dx\
   -(gamma*(u1*u1*u2-u1+a))*q*dx \
  + ((u2 - u20) / k)*v*dx + d2*inner(grad(u2), grad(v))*dx\
  -(gamma*(-u1*u1*u2+b))*v*dx 

ax = dot((Xh / k),n_oriented)*phi*dx  \
     + dot(H*n_oriented,Xi)*dx - inner(grad(Xh),grad(Xi))*dx

Lx = dot(kp,u_vec)*phi*dx + dot((Xh0 / k),n_oriented)*phi*dx

# Create nonlinear problem and Newton solver

problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver

ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

#Linear problem

problemX = LinearProblem(ax, Lx)

# Step in time
t = 0.0

# Output file
name = "mesh_at_t"+str(t)
mesh.name = name

# Get the sub-space for u1 and X and the corresponding dofs 

V0, dofs = ME.sub(0).collapse()

u1 = u.sub(0)
u2 = u.sub(1)

file1 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/outputu1.xdmf", "w")
file2 = XDMFFile(MPI.COMM_WORLD, "SchnakenbergProof/H.xdmf", "w")

file1.write_mesh(mesh)
file2.write_mesh(mesh)

file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file2.write_function(H0, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

norm_stop = 0.1e-6

fileL2 = open('L2.dat', 'w') 

z = 0
inc = 0 
impresion = 1

u0.x.array[:] = u.x.array

t += dt

X_h = problemX.solve()
r = solver.solve(u)

u0.x.array[:] = u.x.array

Xn, Hn = X_h.sub(0).collapse(), X_h.sub(1).collapse()

Xh0.x.array[:] = Xn.x.array[:]
H0.x.array[:] = Hn.x.array[:]

X_coord2 = Xh0.x.array.reshape(len(mesh.geometry.x[:,1]),2)
xf2 = Xh0.x.array[:]

mesh.geometry.x[:,0] = X_coord2[:,0]
mesh.geometry.x[:,1] = X_coord2[:,1]

# Save solution to file (VTK)

name = "mesh_at_t"+str(t)
mesh.name = name

file1.write_mesh(mesh)
file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file2.write_mesh(mesh)
file2.write_function(Xh0, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

z += 1

file1.close()
file2.close()

fileL2.close()

Thank you in advance.