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: