Initial Conditions for mixed Element

Martin1 · December 2, 2019, 6:45pm

Hello Everyone,

I would like to define Inital Conditions for a mixed Element Problem. I have found this post: https://fenicsproject.org/qa/13118/specifying-initial-values-for-mixed-vector-and-scalar-fields/
and tried to implement it in my Code, which runs without errors. But when I compare the results with the ones, were I dont have the Inital Conditions they are the same. (I used tkdiff and there was no singel diffrence in the files I looked at).
I would really appreciate it, if you could tell me what I did wrong.
Here is a running example: I hope its not to long

from fenics import *
import matplotlib.pyplot as plt

Elasticity parameters

E= 10010**6 # E in Pa
nu = 0.25 # Querkontraktionszahl
n_por = 0.01 # Porosität wird zu 0 gesetzt da sie bei ihnen garnicht erst vorkommt
k0= 110**(-14)
visk = 1 * 10**(-5)
k = k0/visk # Permeabilität in m^4/Ns
b = 1 # Biot Koefficent
M = 10**14 # Biot Modul
F = 1000 # Kraft welche auf die Platte wirkt
a = 1 # Länge des Plattenabschnitts
G = Constant(E/(2*(1 + nu))) # Schubmodul
K = Constant((1+nu)/(3*(1-nu))2G*(1-nu)/(1-2nu))
B_big = Constant(bM/(K+pow(b,2)M))
nu_u = Constant((3nu+bB_big(1-2nu))/(3-bB_big*(1-2*nu)))

Zeitvariablen

T = 1.0 # final time
num_steps = 1000 # number of time steps
dt = T / num_steps # time step size
#dt = 0.01

Create mesh and define function space

nx0 = 0
ny0 = 0
nx1 = 1
ny1 = 1

mesh = RectangleMesh(Point(nx0, ny0), Point(nx1, ny1), 30,30)
P2 = VectorElement(‘Lagrange’, mesh.ufl_cell(), 2) # für die Verschiebung quadratisch
P1 = FiniteElement(‘Lagrange’, mesh.ufl_cell(), 1) # für den Druck linear

Make a mixed space V.sub(0)- FiniteElement

V.sub(1) - VectorElement

TH = P2 * P1
V = FunctionSpace(mesh, TH)

Mark boundary subdomians

class rechts(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],1,tol)
class links(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],0,tol)
class oben(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1],1,tol)
class unten(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1],0,tol)

right = rechts()
left = links()
top = oben()
bottom = unten()

boundaries = MeshFunction(‘size_t’,mesh,mesh.topology().dim()-1)
boundaries.set_all(0)

left.mark(boundaries,1)
right.mark(boundaries,2)
top.mark(boundaries,3)
bottom.mark(boundaries,4)

Define functions

(u, p) = TrialFunctions(V)
(v, r) = TestFunctions(V)
du = Function(V)
u0, p0 = split(du)

Define Constants

n = FacetNormal(mesh)
n_por = Constant(n_por) # Porosität
k = Constant(k) # Permeabilität
nu = Constant(nu) # Querkontraktionszahl
b = Constant(b)
M = Constant(M)
F = Constant(F)
a = Constant(a)
mu = Constant(E/(2*(1 + nu)))
lam = Constant(E * nu/( (1+nu)(1-2 nu)))

Fluidfluss am Rand

q = 0 # -q_vektor*n= q

Kinematics

d = u.geometric_dimension()
I = Identity(d) # Identity tensor

linearer Spannungstensor

def epsilon(u):

return 0.5*(nabla_grad(u) + nabla_grad(u).T)
return sym(nabla_grad(u))
Spur des Spannungstensors

def tr_epsilon(u):
return tr(epsilon(u))

lineare Festkörperspannung

def sigma_s(u):
return lam* tr_epsilon(u)* Identity(d) + 2* mu * epsilon(u)

Spannungstensor aus Festkörper und Fluidanteil

def sigma(u,p):
return (1-n_por) * sigma_s(u) -n_por* p*I

ds_test = ds(domain=mesh, subdomain_data=boundaries)

Impulsbilanz + Darcy und Fluidmassenbilanz

T1 = Expression((“1000”, “0”),degree =2)
T2 = Expression((“-1000”, “0”),degree =2)

Imp =inner(sigma(u,p),nabla_grad(v)) * dx -inner(T1,v) * ds_test(3) -inner(T2,v) * ds_test(4)
Fl = qrds_test(3) +qrds_test(4) - k * inner(nabla_grad( p ),nabla_grad( r )) * dx -
(tr_epsilon(u) - tr_epsilon(u0))/dt *r *dx

eqn = Imp + Fl

Fluid +Verschiebungsrandbedingung - free-slip BC

bcu_B = DirichletBC(V.sub(0).sub(1), Constant((0)), bottom)
bcu_L = DirichletBC(V.sub(0).sub(0), Constant((0)), left)
bcp_R = DirichletBC(V.sub(1), Constant((0)), right) >
bc_ges = [bcu_L, bcu_B, bcp_R]

Initial Conditions:

e_p0 = Expression(‘1/(3a) * B_big * (1+nu_u) * F’, degree =2, F=F, a=a, B_big=B_big,nu_u=nu_u)
e_u0 = Expression(('(Fnu_u/(2*G))x[0]/a’ , '-F(1-2 * nu_u)/(2 * G) * (x[1]/a)’), degree =2, F=F, a=a, G=G ,nu_u=nu_u)

p_0 = interpolate(e_p0, V.sub(1).collapse())
u_0 = interpolate(e_u0, V.sub(0).collapse())
w = Function(V)
assign(w, [u_0, p_0])

Aufteilen der Gleichung in linke und rechte Seite

l = lhs(eqn)
r = rhs(eqn)

Solve variational problem

Step in time

t = 0

for n in range(num_steps):
# Update current time
t += dt

solve(l == r, w, bc_ges)
(u,p) = w.split(True)

assign(du.sub(0), u)

kamensky · December 2, 2019, 7:24pm

If I’m interpreting the variables correctly, the previous time step’s solution is du, so that is what should be set to the initial condition, not w, which does not appear in the problem eqn (and therefore does not affect its solution), and is simply over-written with the current solution after the solve.

Martin1 · December 2, 2019, 8:05pm

Thank you so much for your quick and good answer!
I changed
assign(w, [u_0, p_0])
to
assign(du, [u_0, p_0])
and now its working!

Martin1 · February 16, 2020, 9:36am

Hello,

I’m trying to solve a different problem: I have a poromechanic Material and I’m pulling on the right side.

I’m not sure, if my implementation of the initial Condition is correct, because I define, that the pressure at
p0 = 10000
but the pressure of my next time step is in a range from 0 (because of my boundary Condition) to -100. I dont understand the big difference between t0 and t1. Is there a mistake in my implementation? Or do I have to implement an initial pressure differently?

Here is the section with the initial condition of my code (and at the bottom is a running example)

P2 = VectorElement('Lagrange', Netz.ufl_cell(), 2)      # displacement
P1 = FiniteElement('Lagrange', Netz.ufl_cell(), 1)      # pressure

# Make a mixed space 
TH = P2 * P1
V = FunctionSpace(Netz, TH)


# Class representing the intial conditions
class InitialConditions(UserExpression):
    def eval(self, values, x):
        values[0] = 0.0                 # ux
        values[1] = 0.0                 # uy
        values[2] = 10000.0             # p
    def value_shape(self):
        return (3,)
# Define functions

(u, p) = TrialFunctions(V)		# new Timestep
(v, r)  = TestFunctions(V)           
du = Function(V)                # old Timestep
u0, p0 = split(du)

# Create intial conditions and interpolate
u_init = InitialConditions()
du.interpolate(u_init)










'''2D Poromechanik'''
from fenics import *
from numpy import cos , pi
from ufl import nabla_div


parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True

# Elasticity parameters
E= 220*1000		# E in Pa
nu = 0.45 		# Querkontraktionszahl
b = 1                # Biot Koeffizient b = 1 für inkompr. Fluid + Solid (Castelletto)
k = 2.88*10**(-14)		# Permeabilität	in m^4/Ns
K_f = 3.3*10**9 #von Wasser in Pa Detournay S- 24
phi = 0.54 # Porosität
M = Constant(K_f/phi)
p_ini = 10000

# Timevariables
T = 10          # final time
num_steps = 100    # number of time steps
dt = T / num_steps # time step size

# Define Constants
b = Constant(b)		# Biot Koeff
k   = Constant(k)			# Permeabilität 
nu  = Constant(nu)			# Querkontraktionszahl
G = Constant(E/(2*(1 + nu)))
lam = Constant(E*nu/( (1+nu)*(1-2*nu)))
K = Constant((1+nu)/(3*(1-nu))*2*G*(1-nu)/(1-2*nu))
B_big = Constant(b*M/(K+pow(b,2)*M))
nu_u = Constant((3*nu+b*B_big*(1-2*nu))/(3-b*B_big*(1-2*nu)))



##############################################################################
#                                   Create mesh  Beginning                      #
##############################################################################

# Create mesh and define function space 
nx0 = 0
ny0 = 0
nx1 = 0.090 # in m
ny1 = 0.030 # in m
#nx1 = 9*10^-3
#ny1 = 3*10^-3

def mesh(lx,ly, Nx,Ny):
    m = UnitSquareMesh(Nx, Ny)
    x = m.coordinates()

    #Refine on the left and right sides
    x[:,0] = (x[:,0] - 0.5) * 2.
    x[:,0] = 0.5 * (cos(pi * (x[:,0] - 1.) / 2.) + 1.)   

    #Scale
    x[:,0] = x[:,0]*lx
    x[:,1] = x[:,1]*ly

    return m

Netz = mesh(nx1,ny1, 40,40)
nCells = Netz.num_cells()
print('Das Netz besteht aus', nCells,'Elementen')

##############################################################################
#                                   Create mesh  Ende                      #
##############################################################################



P2 = VectorElement('Lagrange', Netz.ufl_cell(), 2)      # displacement
P1 = FiniteElement('Lagrange', Netz.ufl_cell(), 1)      # pressure

# Make a mixed space 
TH = P2 * P1
V = FunctionSpace(Netz, TH)

##############################################################################
#                                Boundary Stuff: Anfang                      #
##############################################################################
# Mark boundary subdomians
def boundary(x, on_boundary):   # fuer Druck-RB
    return on_boundary
# Mark boundary subdomians
class rechts(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],nx1,tol)
class links(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],nx0,tol)
class oben(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1],ny1,tol)
class unten(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1],ny0,tol)
    
right = rechts()
left = links()
top = oben()
bottom = unten()

boundaries = MeshFunction('size_t',Netz,Netz.topology().dim()-1)

boundaries.set_all(0)

left.mark(boundaries,1)
right.mark(boundaries,2)
top.mark(boundaries,3)
bottom.mark(boundaries,4)

ds_test = ds(domain=Netz, subdomain_data=boundaries)
n   = FacetNormal(Netz)

# Boundary Conditions
bcp_O  = DirichletBC(V.sub(1), Constant(0), top)   
bcp_U  = DirichletBC(V.sub(1), Constant(0), bottom)   

bcu_L  = DirichletBC(V.sub(0), Constant((0, 0)), left)

# pulling on the right side
Dehnrate = 1*10**(-3) 
u_RB_lin = Expression(("Dehnrate*t"), degree=2, Dehnrate=Dehnrate, t=0)
RB = DirichletBC(V.sub(0).sub(0),  u_RB_lin, right)

bc_ges = [bcu_L, bcp_U, bcp_O, RB]

# Class representing the intial conditions
class InitialConditions(UserExpression):
    def eval(self, values, x):
        values[0] = 0.0                 # ux
        values[1] = 0.0                 # uy
        values[2] = 10000.0             # p
    def value_shape(self):
        return (3,)
##############################################################################
#                                Boundary Stuff: Ende                        #
##############################################################################
# Define functions
Tensor_grad = 2

(u, p) = TrialFunctions(V)		# new Timestep
(v, r)  = TestFunctions(V)           
du = Function(V)                # old Timestep
u0, p0 = split(du)

# Create intial conditions and interpolate
u_init = InitialConditions()
du.interpolate(u_init)


theta = 0.5     # 1 -> Euler Rückwärts, 0.5 -> Crank Nicolson
p_mid = (1-theta)*p0+ theta*p

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor

# linearer strain tensor
def epsilon(u):
#    return 0.5*(nabla_grad(u) + nabla_grad(u).T) 
    return sym(nabla_grad(u))

def tr_epsilon(u):
    return tr(epsilon(u))    

# Stress tensor
def sigma(u,p):
	return 2*G* epsilon(u) + 2*G * nu/(1-2*nu) *(nabla_div(u))*I - b*p*I

# Momentum balance + Darcy und Fluidmassbalance	
Imp =inner(sigma(u,p),nabla_grad(v))*dx 
Fl = b*(tr_epsilon(u) - tr_epsilon(u0))/dt *r *dx + 1/M * (p-p0)/dt *r*dx +\
     k*inner(nabla_grad((p_mid)),nabla_grad(r))*dx 

F = Imp + Fl

l = lhs(F)
r = rhs(F)

t = 0

##############################################################################
#                               Lösungsschleife                        #
##############################################################################       
for n in range(num_steps):
    # Update current time
    t += dt
    u_RB_lin.t=t


    w = Function(V)
    if n ==0:
        w.interpolate(u_init)
 
    solve(l == r, w, bc_ges)
    (u,p) = w.split(True)
    du.assign(w)

dokken · February 16, 2020, 10:03am

as far as I can tell, you are not updating du (u0, p0) with the previous solution in your loop. You should do du.assign(w)

Martin1 · February 16, 2020, 10:48am

ah, Iam sorry I deleted the whole saving of data lines after the splitting to make the code shorter and also deleted the line:

du.assign(w)

but it is in my code, so unfortunatly, its not that, but thank you for your quick answer.
(I will edit my previous post to put it in there)

Martin1 · February 17, 2020, 3:23pm

Can I also assign a value to a TrialFunction, like I do it for a “normal” function like I do it here:
du.interpolate(u_init) ?

Because something like p.assign(…) failed.

dokken · February 17, 2020, 7:01pm

A trial function can’t be assigned values

Kart · August 12, 2020, 10:30am

Hi, I have 2 questions.
For example in my problem:

P3 = FiniteElement('P', mesh.ufl_cell(), 2)
ele = MixedElement([P3,P3,P3])
VV = FunctionSpace(mesh, ele)
#Define test functions
vv1,vv2,vv3 = TestFunction(VV)
#Define problem to solve
uu = Function(VV)
phi,bound_s,bound_ns = split(uu)
#Define initial and boundary conditions
u_n = Function(VV) #Automatically sets to zero
phi_n,bound_s_n,bound_ns_n = split(u_n)

1) How do I assign an inital value to phi_n?

Given that I have the defined various boundaries and regions within my mesh as follows:

mesh = Mesh()
with XDMFFile(output_folder_name+model+".xdmf") as infile:
    infile.read(mesh) #Feeds data from xdmf to mesh variable

mvc = MeshValueCollection("size_t", mesh, 1) #Container for mesh values
with XDMFFile(output_folder_name+"mf"+model+".xdmf") as infile:
    infile.read(mvc, "name_to_read") #Feeds boundary data to mvc 

mvc2 = MeshValueCollection("size_t", mesh, 1) 
with XDMFFile(output_folder_name+"cf"+model+".xdmf") as infile:
    infile.read(mvc2, "name_to_read") 

boundaries = MeshFunction("size_t", mesh, mvc) 
regions = MeshFunction("size_t", mesh, mvc2)

2) How do I assign an inital value to phi_n only for a specific region or boundary?

Please advise, thanks

dokken · August 12, 2020, 10:45am

For the first question, see: Access function in time dependent XDMFFile which will refer you to the FunctionAssigner-function.
For the second question, consider this minimal example:

from dolfin import *
mesh = UnitSquareMesh(2,1)
mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
class left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]<1e-10 and on_boundary

left().mark(mf, 1)

P3 = FiniteElement('P', mesh.ufl_cell(), 2)
ele = MixedElement([P3,P3,P3])
VV = FunctionSpace(mesh, ele)
#Define test functions
vv1,vv2,vv3 = TestFunction(VV)
#Define problem to solve
uu = Function(VV)
phi,bound_s,bound_ns = split(uu)
#Define initial and boundary conditions
u_n = Function(VV) #Automatically sets to zero
phi_n,bound_s_n,bound_ns_n = split(u_n)


V0,V1,V2 = VV.split()
V0 = V0.collapse()
v0 = Function(V0)
v0.interpolate(Expression("1-x[0]+x[1]",degree=1))

# Setting values for part of boundy of subspace_vector
init_bc = DirichletBC(VV.sub(0), v0, mf, 1)
init_bc.apply(u_n.vector())
print("u_n", u_n.vector().get_local())
print("sub-vector:", u_n.split(deepcopy=True)[0].vector().get_local())
# Setting bc for whole boundary
full_bc = DirichletBC(VV, Expression(("1-x[0]","1-x[0]","1-x[0]"),degree=1), "on_boundary")
full_bc.apply(u_n.vector())
print("Full space", u_n.vector().get_local())

Kart · August 12, 2020, 11:50am

Hi thanks for getting back to me. Sorry but I don’t think I was clear with my question. I don’t mean to apply a Dirichlet boundary condition but rather an initial condition.
Essentially, I want to assign a value to phi_n for a specific edge which is defined in boundaries at time = 0.
Where
phi_n, bound_s_n, bound_ns_n = split(u_n)
and
boundaries = MeshFunction("size_t", mesh, mvc)
such that for time > 0 the values at this edge changes according to the diffusivity of the material.

dokken · August 12, 2020, 11:52am

That is what i have done here. I am just using the Dirichlet BC to apply the specific condition to the vector at the given edges.
After applying the BC to u_n, you can use it as an initial condition in the variational form.

Kart · August 12, 2020, 11:55am

Oh I see. Thank you very much!

Topic		Replies	Views
How to define initial conditions in mixed formulation using dolfinx? dolfinx	3	1082	June 24, 2022
Problem defining initial conditions : Expression problem? I/O	2	3942	April 18, 2019
How to set the initial conditions in FEniCSx General	2	936	December 7, 2022
Using the solution of Stokes (mixed elements) as initial condition for Chorin (not mixed)	3	725	January 22, 2021
Defining initial condition on 1-D IntervalMesh I/O	1	625	February 2, 2022

Initial Conditions for mixed Element

Elasticity parameters

Zeitvariablen

Create mesh and define function space

Make a mixed space V.sub(0)- FiniteElement

V.sub(1) - VectorElement

Mark boundary subdomians

Define functions

Define Constants

Fluidfluss am Rand

Kinematics

linearer Spannungstensor

return 0.5*(nabla_grad(u) + nabla_grad(u).T)

Spur des Spannungstensors

lineare Festkörperspannung

Spannungstensor aus Festkörper und Fluidanteil

Impulsbilanz + Darcy und Fluidmassenbilanz

Fluid +Verschiebungsrandbedingung - free-slip BC

Initial Conditions:

Aufteilen der Gleichung in linke und rechte Seite

Solve variational problem

Step in time

Related Topics