Using the solution of Stokes (mixed elements) as initial condition for Chorin (not mixed)

1

Hi,

I get an error if I try to import and use as initial condition my velocity field obtained by solving a mixed problem (the Stokes Problem). I want to use this initial condition to solve the non mixed problem of Navier.Stokes through a projection mehtod (non mixed space). How can I solve this?
Here the error I get an the code:

"from future import print_function
from fenics import *
from ufl import nabla_grad
from ufl import nabla_div
import matplotlib.pyplot as plt
from mshr import *
import numpy as np

mesh = Mesh(“cylinder.xml.gz”)
n = FacetNormal(mesh)

#Create space
P1 = FiniteElement(‘Lagrange’, triangle, 1)
P2 = VectorElement(‘Lagrange’,triangle,2)
element = MixedElement([P2,P1])
W = FunctionSpace(mesh, element)

w = Function(W)

#IMPORT SOLUTION
timeseries_w = TimeSeries(“velocity_series”)
timeseries_w.retrieve(w.vector(), 1)
u_n, p_n = w.split()

#SET PARAMETERS
T = 1.5 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
k = Constant(dt)
vu = 1 # kinematic viscosity mu/rho
vu = Constant(vu)
f = Constant((0, 0))

#Set Fuction for Chorin
V = VectorFunctionSpace(mesh, ‘P’, 2)
Q = FunctionSpace(mesh, ‘P’, 1)

Define trial and test functions

u_t = TrialFunction(V)
v = TestFunction(V)
p_t = TrialFunction(Q)
q = TestFunction(Q)

Define functions for solutions

u_ = Function(V)
p_ = Function(Q)

#Define boundaries
up_down = ‘near(x[1], 1)||near(x[1], -1)’
links = ‘near(x[0], -1)’
rechts=‘near(x[0], 4)’
zero_con=‘on_boundary && x[0]>-0.3 && x[0]<0.3 && x[1]>-0.3 && x[1]<0.3’

And b.c.

f2 = Expression((“1-x[1]*x[1]”, “0.0”), degree=2)
g=0.5 #Neumann non homogeneous
g = Constant(g)
bc_up_down = DirichletBC(V, Constant((0.0, 0.0)), up_down)
bc_kreis= DirichletBC(V, Constant((0.0, 0.0)), zero_con)
bc_left=DirichletBC(V, f2, links)
bcu=[bc_up_down,bc_kreis,bc_left]

Define strain-rate tensor

def epsilon(u):
return sym(nabla_grad(u))

Define stress tensor

def sigma(u, p):
return 2vuepsilon(u) - pIdentity(len(u)) #p/rho zu ergänzen (auch in Randterme)
U = 0.5(u_n + u_t)

#Define Variational problem
F1 = dot((u_t - u_n) / k, v)*dx +
dot(dot(u_n, nabla_grad(u_n)), v)*dx \

  • inner(sigma(U, p_n), epsilon(v))*dx \
  • dot(p_n*n, v)ds - dot(vunabla_grad(U)*n, v)ds-dot(gn,v) \
  • dot(f, v)*dx
    a1 = lhs(F1)
    L1 = rhs(F1)"

I then get the error:

grafik
grafik1044×456 54.7 KB

2

Please follow the guidelines Read before posting: How do I get my question answered?
and format all code by encapsulating it with ```.
Also remove all code node required to reproduce the error. It would be Good if you could reproduce the problem with the built in meshes.
As you are trying to take a function from a Mixed space and assign it to a non Mixed problem, you can use
u,p=w.split(deepcopy=True)

1 Like
3

Thank you Dokken the extended split command seems to be what I searched for. I tried to upload the gmz and the timeseries, but it was not possible (unauthorized). Hence I can quote here form the generatin code:
``# Create mesh and normals
L=4
W=1
rechteck = Rectangle(Point(-1, -1), Point(L, W))#, 10, 10, "right/left")
kreis = Circle(Point(0, 0), 0.2)
mesh_form = rechteck- kreis
mesh=generate_mesh(mesh_form,20)
#Create space
P1 = FiniteElement('Lagrange', triangle, 1)
P2 = VectorElement('Lagrange',triangle,2)
element = MixedElement([P2,P1])
V = FunctionSpace(mesh, element)
u__ = Function(V)

I then solve Stokes for a u__
solve(a == L, u__, bcu)
File("cylinder.xml.gz") << mesh
timeseries_u = TimeSeries("velocity_series")
timeseries_u.store(u__.vector(),1)

4

Please format all your code using ```

Such that it renders properly:

from dolfin import *
mesh=UnitSquareMesh(5,5)
...