Giving initial condition in FEniCSx0.9.0 using mixed_element

dolfinx
1

Hi, everyone! I’m tring to solving the swelling problem of gels using the mixed space. My codes are as follows and I found that The three-directional deformation gradient component F33 in the problems of displacement, chemical potential and plane stress cannot be correctly initialized：

# Function space

V_dis = element("CG",domain.basix_cell(), 2,shape=(2,)) #displacement
V_mu = element("CG",domain.basix_cell(), 1)          #chemical potential   
V_F33 = element("CG",domain.basix_cell(), 1)     #P33 
V_Mix = mixed_element([V_dis, V_mu, V_F33])
Mix = dfx.fem.functionspace(domain, V_Mix)

dis_space = Mix.sub(0)  # 位移空间(向量)
V_dis, _ = dis_space.collapse()
mu_space = Mix.sub(1)   # 化学势空间
V_mu, _ = mu_space.collapse()
F33_space = Mix.sub(2)  # F33空间
V_F33, _ = F33_space.collapse()

# Functions
u = dfx.fem.Function(Mix)        # dis,mu,F33 in current time step
dis, mu, F33 = split(u)

u_old = dfx.fem.Function(Mix)    # dis,mu,F33 in last time step
dis_old, mu_old, F33_old = split(u_old)

u_test = ufl.TestFunction(Mix)
dis_test, mu_test, F33_test = split(u_test)

du = ufl.TrialFunction(Mix)
ddis, dmu, dF33 = split(du)

u.sub(0).x.array[:] = 0.0

u.sub(1).x.array[:] = muini.value

u.sub(2).x.array[:] = 1.0

u_old.sub(0).x.array[:] = 0.0

u_old.sub(1).x.array[:] = muini.value

u_old.sub(2).x.array[:] = 1.0

print(f"Initial u.sub(0): {u.sub(0).x.array}")

print(f"Initial u.sub(1): {u.sub(1).x.array}")

print(f"Initial u.sub(2): {u.sub(2).x.array}")

print(f"Initial u_old.sub(0): {u_old.sub(0).x.array}")

print(f"Initial u_old.sub(1): {u_old.sub(1).x.array}")

print(f"Initial u_old.sub(2): {u_old.sub(2).x.array}")

and the output initial value is：
Initial u.sub(0): [1.84369148e-05 1.84369148e-05 1.84369148e-05 … 1.84369148e-05
1.84369148e-05 1.84369148e-05]
Initial u.sub(1): [1.84369148e-05 1.84369148e-05 1.84369148e-05 … 1.84369148e-05
1.84369148e-05 1.84369148e-05]
Initial u.sub(2): [1.84369148e-05 1.84369148e-05 1.84369148e-05 … 1.84369148e-05
1.84369148e-05 1.84369148e-05]
Initial u_old.sub(0): [1. 1. 1. … 1. 1. 1.]
Initial u_old.sub(1): [1. 1. 1. … 1. 1. 1.]
Initial u_old.sub(2): [1. 1. 1. … 1. 1. 1.]
can anyone help me, thanks!

2

the sub(i) on a mixed space only gives a window into that entry. The x vectors still map to the original full space. To access only subspaces of that vector, you need precisely the return argument that you crossed out in the collapse action:

V_dis, V0_map = dis_space.collapse()
u.x.array[V0_map] = ...

So the V0_map array.

P.s. please make code reproducible (including includes, mesh, etc), so its easy to copy/paste. That helps a lot with increasing the efficiency and quality of responses.

3

Thanks your reply! It seems that this operation dose not work. The initial values are just applied on the last three dofs:
‘’'# Function space

V_dis = element(“CG”,domain.basix_cell(), 2,shape=(2,)) #displacement
V_mu = element(“CG”,domain.basix_cell(), 1) #chemical potential
V_F33 = element(“CG”,domain.basix_cell(), 1) #P33
V_Mix = mixed_element([V_dis, V_mu, V_F33])
Mix = dfx.fem.functionspace(domain, V_Mix)

dis_space = Mix.sub(0) # 位移空间(向量)
V_dis, _ = dis_space.collapse()
mu_space = Mix.sub(1) # 化学势空间
V_mu, _ = mu_space.collapse()
F33_space = Mix.sub(2) # F33空间
V_F33, _ = F33_space.collapse()

Functions

u = dfx.fem.Function(Mix) # dis,mu,F33 in current time step
dis, mu, F33 = split(u)

u_old = dfx.fem.Function(Mix) # dis,mu,F33 in last time step
dis_old, mu_old, F33_old = split(u_old)

u_test = ufl.TestFunction(Mix)
dis_test, mu_test, F33_test = split(u_test)

du = ufl.TrialFunction(Mix)
ddis, dmu, dF33 = split(du)

dis_col, dis_col_to_Mix = Mix.sub(0).collapse()

u_old.x.array[dis_col_to_Mix] = 0.0

u.x.array[dis_col_to_Mix] = 0.0

u_old.x.scatter_forward ()

# chemical potential

mu_col, mu_col_to_Mix = Mix.sub(1).collapse()

u_old.x.array[mu_col_to_Mix] = muini.value

u.x.array[mu_col_to_Mix] = muini.value

u_old.x.scatter_forward ()

# F33

F33_col, F33_col_to_Mix = Mix.sub(2).collapse()

u_old.x.array[F33_col_to_Mix] = 1.0

u.x.array[F33_col_to_Mix] = 1.0

u_old.x.scatter_forward ()‘’’
and the output are:
‘’’
Initial u.sub(0): [0.00000000e+00 0.00000000e+00 0.00000000e+00 … 0.00000000e+00
1.84369148e-05 1.00000000e+00]
Initial u.sub(1): [0.00000000e+00 0.00000000e+00 0.00000000e+00 … 0.00000000e+00
1.84369148e-05 1.00000000e+00]
Initial u.sub(2): [0.00000000e+00 0.00000000e+00 0.00000000e+00 … 0.00000000e+00
1.84369148e-05 1.00000000e+00]
Initial u_old.sub(0): [0.00000000e+00 0.00000000e+00 0.00000000e+00 … 0.00000000e+00
1.84369148e-05 1.00000000e+00]
Initial u_old.sub(1): [0.00000000e+00 0.00000000e+00 0.00000000e+00 … 0.00000000e+00
1.84369148e-05 1.00000000e+00]
Initial u_old.sub(2): [0.00000000e+00 0.00000000e+00 0.00000000e+00 … 0.00000000e+00
1.84369148e-05 1.00000000e+00]
‘’’

4

As already mentioned, making your example reproducible makes it easier for us to give you guidance. Please also use 3x` encapsulation for code to make it readable and formatted correctly.

5

Here are my codes,thanks!

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx as dfx
import dolfinx.fem.petsc as petsc
from dolfinx.fem import Constant
from dolfinx.nls.petsc import NewtonSolver
import ufl
from basix.ufl import element, mixed_element
import time
from ufl import ln, inner, grad, tr, Identity, det, dot, inv, exp, div, diff, variable, sinh, derivative, split
from dolfinx.io import XDMFFile
from datetime import datetime

# Set level of detail for log messages (integer)
comm = MPI.COMM_WORLD
mpiRank = MPI.COMM_WORLD.rank
# log.set_log_level(log.LogLevel.WARNING)


##################################################################################
###                             Define Geometr Mesh                            ###
##################################################################################
L0 = 0.01  # length
H0 = 0.01  # height
domain = dfx.mesh.create_rectangle(comm, [[0.0, 0.0], [L0, H0]], [20, 20],
                                   cell_type=dfx.mesh.CellType.quadrilateral)
x = ufl.SpatialCoordinate(domain)

topology_dim = domain.topology.dim

domain.topology.create_connectivity(topology_dim-1, topology_dim)

##################################################################################
###            Functions for Finding Different Areas, Boundary              ###
##################################################################################
#                    4
#         ----------------------
#       1 |                    | 3
#         |                    |
#         ----------------------
#                    2
def Left(x):
    return np.isclose(x[0], 0.0)

def Right(x):
    return np.isclose(x[0], L0)

def Bottom(x):
    return np.isclose(x[1], 0.0)

def Top(x):
    return np.isclose(x[1], H0)

boundaries = [(1, Left),(2,Bottom),(3,Right),(4,Top)]

facet_boundary, facet_markers = [], []        # boundary index 
fdim = topology_dim-1                         # Facet dim
for (marker, locator) in boundaries:
    facets = dfx.mesh.locate_entities(domain, fdim, locator)
    facet_boundary.append(facets)
    facet_markers.append(np.full_like(facets, marker))

facet_boundary = np.hstack(facet_boundary).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_boundary)

facet_tag = dfx.mesh.meshtags(domain, fdim, facet_boundary[sorted_facets], facet_markers[sorted_facets])


##################################################################################
###                   Simulation time-control releated parameters              ###
##################################################################################
# Loading parameters
t = 0.0                 # Initialization of time
dt = 1                 # Fixed step size
numSteps = 4000         # Number of steps
Ttol = numSteps*dt      # Total simulation time
cnt = 0                 # iters
save_cnt = 1            # save_cnt iters save one result
vload = 1e-4            # max dfdsisplacement BC


##################################################################################
###                              Materials Parameters                          ###
##################################################################################
# Material Constants
Nv0 = Constant(domain,PETSc.ScalarType(2.6e-4))    # density of chains
v = Constant(domain,PETSc.ScalarType(1e-28))       # volume of per water molecules
N0 = Constant(domain,PETSc.ScalarType(Nv0/v))
kT = Constant(domain,PETSc.ScalarType(4.05e-21))   # Boltzmann's constant
chi = Constant(domain,PETSc.ScalarType(0.481))     # Flory-Huggins parameter
Dc = Constant(domain,PETSc.ScalarType(1.0e-7))     # Diffusivity
b = Constant(domain,PETSc.ScalarType(2e-9))        # Length of Kuhn segments  
lamda0 = Constant(domain,PETSc.ScalarType(2.78))   # Initial swelling stretch
n = Constant(domain,PETSc.ScalarType(542))         # Number of Kuhn segments

beita0=Constant(domain,PETSc.ScalarType(lamda0*1*(n**0.5)/n*(3-(lamda0**2)*(1**2)*n/(n**2))/(1-(lamda0**2)*(1**2)*n/(n**2))))

muini = Constant(domain,PETSc.ScalarType(Nv0*(n**0.5)/(3.*lamda0**2.)*beita0+ln(1-1/lamda0**3.)+1/lamda0**3.+chi/lamda0**6.))     # mu0/kT

print('beita0:',beita0.value,'mu0/kT:', muini.value)


##################################################################################
###                    Define Function space and Functions                     ###
##################################################################################

# Function space

V_dis = element("CG",domain.basix_cell(), 2,shape=(2,)) #displacement
V_mu = element("CG",domain.basix_cell(), 1)          #chemical potential   
V_F33 = element("CG",domain.basix_cell(), 1)     #P33 
V_Mix = mixed_element([V_dis, V_mu, V_F33])
Mix = dfx.fem.functionspace(domain, V_Mix)

dis_space = Mix.sub(0)  # 位移空间(向量)
V_dis, _ = dis_space.collapse()
mu_space = Mix.sub(1)   # 化学势空间
V_mu, _ = mu_space.collapse()
F33_space = Mix.sub(2)  # F33空间
V_F33, _ = F33_space.collapse()

# Functions
u = dfx.fem.Function(Mix)        # dis,mu,F33 in current time step
dis, mu, F33 = split(u)

u_old = dfx.fem.Function(Mix)    # dis,mu,F33 in last time step
dis_old, mu_old, F33_old = split(u_old)

u_test = ufl.TestFunction(Mix)
dis_test, mu_test, F33_test = split(u_test)

du = ufl.TrialFunction(Mix)
ddis, dmu, dF33 = split(du)

##################################################################################
###                                Initial conditions                          ###
##################################################################################
# displacement
dis_col, dis_col_to_Mix = Mix.sub(0).collapse()
u_old.x.array[dis_col_to_Mix] = 0.0
u.x.array[dis_col_to_Mix] = 0.0
u_old.x.scatter_forward ()
# # chemical potential
mu_col, mu_col_to_Mix = Mix.sub(1).collapse()
u_old.x.array[mu_col_to_Mix] = muini.value
u.x.array[mu_col_to_Mix] = muini.value
u_old.x.scatter_forward ()
# # F33
F33_col, F33_col_to_Mix = Mix.sub(2).collapse()
u_old.x.array[F33_col_to_Mix] = 1.0
u.x.array[F33_col_to_Mix] = 1.0
u_old.x.scatter_forward ()

print(f"Initial u.sub(0): {u.sub(0).x.array}")
print(f"Initial u.sub(1): {u.sub(1).x.array}")
print(f"Initial u.sub(2): {u.sub(2).x.array}")
print(f"Initial u_old.sub(0): {u_old.sub(0).x.array}")
print(f"Initial u_old.sub(1): {u_old.sub(1).x.array}")
print(f"Initial u_old.sub(2): {u_old.sub(2).x.array}")


##################################################################################
###                Kinematics, Free Energy, Sress Tensor, Weak Form            ###
##################################################################################
def r(t):
    return float((float(muini)-0.1)*exp(-t/numSteps)+0.1)

Time_cons = dfx.fem.Constant(domain, PETSc.ScalarType(r(t)))

# Kinematics
d = len(dis)                       # number of dimensions (d=2 for 2D)
I = Identity(d)                    # identity tensor
F = variable(I + grad(dis))        # deformation gradient
F_old = I + grad(dis_old)
J = variable(det(F))*F33           # Jacob tensor
J_old = det(F_old)*F33_old
C = F.T * F
I1 = tr(C) + F33**2.
lamda_chain = (I1/3.)**0.5
beita = lamda0*lamda_chain*(n**0.5)/n*(3.-(lamda0**2.)*(lamda_chain**2.)*n/(n**2.))/(1.- (lamda0**2.)*(lamda_chain**2.)*n/(n**2.))

# Normalized free energy density Wv/kT
Ws = (1/lamda0**3.)*Nv0*n*(lamda0*lamda_chain/(n**0.5)*beita+ln(beita/sinh(beita)))
Wm = (J-lamda0**(-3.))*ln((J-lamda0**(-3.))/J)+((lamda0**3.)*J-1.)*chi/(lamda0**6.)/J-mu*(J-lamda0**(-3.))
W_total = Ws +Wm

# Normalized stress tensor Pv/kT
# PKs = diff(Ws, F)
# PKm = diff(Wm, F)
# PK = PKs +PKm
PK = Nv0/(lamda0**2.)*(n**0.5)*beita/((3.*I1)**0.5)*F + (inv(F).T)*(J*ln((J-lamda0**(-3.))/J)+1/(lamda0**3.)+chi/(lamda0**6.)/J-mu*J)
PK33 = Nv0/(lamda0**2.)*(n**0.5)*beita/((3.*I1)**0.5)*F33 + 1/F33*(J*ln((J-lamda0**(-3.))/J)+1/(lamda0**3.)+chi/(lamda0**6.)/J-mu*J)



# Weak form 
# LinerProblem: _L(terms already know):constant part, _a(including the term unknow): 
dx = ufl.Measure('dx', domain=domain, metadata={'quadrature_degree': 4})
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag,metadata={'quadrature_degree': 4})
Normal_vector = ufl.FacetNormal(domain)

Weak_PK = inner(PK, grad(dis_test))*dx
Weak_PK33 = PK33*F33_test*dx
# Weak_mu = (lamda0**3.)*(J-J_old)*mu_test*dx + dt*Dc*((lamda0**3.)*J-1.)/(lamda0**2.)*dot(dot(inv(C), grad(mu)), grad(mu_test))*dx
Weak_mu = (lamda0**3.)*(J-J_old)*mu_test*dx + dt*Dc*((lamda0**3.)*J-1.)/(lamda0**2.)*dot(dot(grad(mu_test),inv(C)), grad(mu))*dx

Weak = Weak_PK + Weak_PK33 + Weak_mu

Jacbian = derivative(Weak, u, du)


##################################################################################
###                                  Boundaries Conditions                     ###
##################################################################################
# Define dofs for the boundary condition
dofs_1 = dfx.fem.locate_dofs_topological(dis_space.sub(0), facet_tag.dim, facet_tag.find(1))  # Left--x
dofs_2 = dfx.fem.locate_dofs_topological(dis_space.sub(1), facet_tag.dim, facet_tag.find(2))  # Bot--y
dofs_3 = dfx.fem.locate_dofs_topological(mu_space, facet_tag.dim, facet_tag.find(1))  # Left--mu
dofs_4 = dfx.fem.locate_dofs_topological(mu_space, facet_tag.dim, facet_tag.find(2))  # Bot--mu
dofs_5 = dfx.fem.locate_dofs_topological(mu_space, facet_tag.dim, facet_tag.find(3))  # Right--mu
dofs_6 = dfx.fem.locate_dofs_topological(mu_space, facet_tag.dim, facet_tag.find(4))  # Top--mu

# Boundary conditions for displacement
zero = dfx.fem.Constant(domain,dfx.default_scalar_type(0))

bcdisL = dfx.fem.dirichletbc(zero, dofs_1, dis_space.sub(0))   # Left--x--fixed
bcdisB = dfx.fem.dirichletbc(zero, dofs_2, dis_space.sub(1))   # Bottom--y--fixed

# Boundary conditions for the chemical potential
# bcmuL = dfx.fem.dirichletbc(Time_cons, dofs_3, V)     # Left--muini
# bcmuB = dfx.fem.dirichletbc(Time_cons, dofs_4, V)     # Bottom--muini
bcmuR = dfx.fem.dirichletbc(Time_cons, dofs_5, mu_space)     # Right--muini
bcmuT = dfx.fem.dirichletbc(Time_cons, dofs_6, mu_space)     # Top--muini

bc_displacement = [bcdisB, bcdisL]
bc_chemical = [bcmuR, bcmuT]
bc_u = bc_displacement + bc_chemical

##################################################################################
###                                Functions for Visualization                 ###
##################################################################################
Vis_dis = element("CG",domain.basix_cell(), 1,shape=(2,)) 
Vector_space = dfx.fem.functionspace(domain,Vis_dis) #Vector function space
Vis_mu = element("CG",domain.basix_cell(), 1)
Scalar_space = dfx.fem.functionspace(domain,Vis_mu)   #Scalar function space


dis_vis = dfx.fem.Function(Vector_space)
dis_vis.name = "displacement"
dis_expr = dfx.fem.Expression(dis, Vector_space.element.interpolation_points())

mu_vis = dfx.fem.Function(Scalar_space)
mu_vis.name = "chemical potential"
mu_expr = dfx.fem.Expression(mu, Scalar_space.element.interpolation_points())

J_vis = dfx.fem.Function(Scalar_space)
J_vis.name = "J"
J_expr = dfx.fem.Expression(J, Scalar_space.element.interpolation_points())

C_vis = dfx.fem.Function(Scalar_space)
C_vis.name = "C"
C_expr = dfx.fem.Expression((J-1)/v, Scalar_space.element.interpolation_points())


def InterpAndSave(t, file):
    dis_vis.interpolate(dis_expr)
    mu_vis.interpolate(mu_expr)
    J_vis.interpolate(J_expr)
    C_vis.interpolate(C_expr)

    file.write_function(dis_vis, t)
    file.write_function(mu_vis, t)
    file.write_function(J_vis, t)
    file.write_function(C_vis, t)

# Nonliner problem for displacement
u_problem = petsc.NonlinearProblem(Weak,u,bc_u,Jacbian)

displacement_solver = NewtonSolver(comm, u_problem )
displacement_solver.max_it = 100
displacement_solver.atol = 1e-8
displacement_solver.rtol = 1e-8
displacement_solver.set_linesearch = "basic"
displacement_solver.convergence_criterion = "incremental"
displacement_solver.report = True
ksp = displacement_solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly" # "preonly" works equally well
opts[f"{option_prefix}pc_type"] = "lu" # do not use 'gamg' pre-conditioner
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
# opts[f"{option_prefix}ksp_max_it"] = 30
ksp.setFromOptions()

problemName = "Gel_swell"
# xdmf = XDMFFile(domain.comm, "/container_zqf/Gel-Swell-Fenicsx/" + problemName + ".xdmf", "w")
xdmf = XDMFFile(domain.comm, "/container_zqf/Gel-Swell-Fenicsx/" + problemName + ".xdmf", "w")
xdmf.write_mesh(domain)
xdmf.rewrite_function_mesh = False
xdmf.flush_output = True
xdmf.functions_share_mesh = True


InterpAndSave(t, xdmf)

print("------------------------------------")
print("Simulation Start")
print("------------------------------------")
startTime = datetime.now()
while (t < Ttol):
    t += dt
    if mpiRank == 0:  # only prints on processor 0 (parallel runs)
        print("\n时间: ", t, "循环:",cnt, "边界mu：", r(t))
        # print("Strain: ", vload*t)

    Time_cons.value = r(t)

    # Solve displacement
    displacement_solver.solve(u)
    u.x.scatter_forward()

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

    cnt = cnt + 1
    max_displacement = np.max(np.linalg.norm(u.sub(0).x.array.reshape(-1, 2), axis=1))
    max_chemical_potential = np.max(u.sub(1).x.array)
    print(f"最大位移: {max_displacement:.3e}")
    print(f"最大化学势: {max_chemical_potential:.3e}")

    if cnt % 1 == 0:
        InterpAndSave(t, xdmf)
endTime = datetime.now()
elapseTime = endTime - startTime
print("------------------------------------------")
print("Elapsed real time:  {}".format(elapseTime))
print("------------------------------------------")
xdmf.close()
6

Atop of being reproducible, it’d also be great if your code snippets could be minimal.

This code snippet captures the problem you’re running into, right?

from mpi4py import MPI
import dolfinx as dfx
from basix.ufl import element, mixed_element

domain = dfx.mesh.create_unit_square(MPI.COMM_WORLD,2,2)

# Function space
V_dis = element("CG",domain.basix_cell(), 2,shape=(2,))
V_mu = element("CG",domain.basix_cell(), 1)
V_Mix = mixed_element([V_dis, V_mu])
Mix = dfx.fem.functionspace(domain, V_Mix)

V_dis, _ = Mix.sub(0).collapse()
V_mu, _ = Mix.sub(1).collapse()

# Functions
u = dfx.fem.Function(Mix)        # dis,mu in current time step
u.sub(0).x.array[:] = 0.0
u.sub(1).x.array[:] = 1.0

print(f"Initial u.sub(0): {u.sub(0).x.array}") # [1,1,1,1,1,1,...]
print(f"Initial u.sub(1): {u.sub(1).x.array}") # [1,1,1,1,1,1,...]
print("are the same")

This is the solution I suggested:

from mpi4py import MPI
import dolfinx as dfx
from basix.ufl import element, mixed_element

domain = dfx.mesh.create_unit_square(MPI.COMM_WORLD,2,2)

# Function space
V_dis = element("CG",domain.basix_cell(), 2,shape=(2,))
V_mu = element("CG",domain.basix_cell(), 1)
V_Mix = mixed_element([V_dis, V_mu])
Mix = dfx.fem.functionspace(domain, V_Mix)

# Get maps to subspaces
V_dis, V0_map = Mix.sub(0).collapse()
V_mu, V1_map = Mix.sub(1).collapse()

# Functions
u = dfx.fem.Function(Mix)
u.x.array[V0_map] = 0.0
u.x.array[V1_map] = 1.0

print(f"Initial u: {u.x.array}") #[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0.  0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.  0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 1.]
print("Has ones and zeros")
1 Like