Hi every body
I am amateur in FEniCS and I am trying to solve mixed navier-stokes system:
ρ(ut + u. ∇u) -ηΔu + ∇p = ρg in Ω
∇.u = 0 in Ω
ρt + u . ∇ρ = 0 in Ω
u(x; 0) = u0 in Ω
ρ(x; 0) = ρ0 in Ω
u(x; t) = ud on ∂Ω
I used newton method and the variational form seems to be correct but I encountered this error:
"In instant.recompile: The module did not compile, see 'E:\Users\sabacmp.instant
\error\instant_module_cce1c1e65526f1f1814e9c100641696ba0c6dddc\compile.log'
Traceback (most recent call last):
File "step2-cahn-hilliard.py", line 115, in
solver.solve(problem, u.vector())
File "step2-cahn-hilliard.py", line 49, in F
assemble(self.L, tensor=b)"
here is my cod:
from dolfin import *
mesh = UnitSquare(50,50)
# Define time-dependent pressure boundary condition
p_in = Expression("sin(3.0*t)", t=0.0)
# define function spaces
V = FunctionSpace(mesh, "Lagrange", 1)
Q = FunctionSpace(mesh, "CG", 1)
S = FunctionSpace (mesh, "CG", 1)
ME = MixedFunctionSpace([V, Q, S])
# Define boundary conditions
noslip = DirichletBC(V, 0,"on_boundary")
inflowp = DirichletBC(Q, p_in, "x[1] > 1.0 - DOLFIN_EPS")
outflowp = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS")
inflows = DirichletBC(S, 1, "x[1] > 0.5 - DOLFIN_EPS")
outflows = DirichletBC(S, 0.5, "x[1] < 0.5 - DOLFIN_EPS")
bcu = [noslip]
bcp = [inflowp, outflowp]
bcs = [inflows, outflows]
bc = [bcu, bcp, bcs]
## Class representing the intial conditions
class InitialConditions(Expression):
def __init__(self):
pass
def eval(self, values, x):
values[0] = 0.0
values[1] = 0.0
if x[0] > 0.5:
values[2] = 1
else:
values[2] = 0.5
def value_shape(self):
return (3,)
# Class for interfacing with the Newton solver
class NavierStokesEquation(NonlinearProblem):
def __init__(self, a, L, bc):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bc = bc
self.reset_sparsity = True
def F(self, b, x):
assemble(self.L, tensor=b)
self.bc.apply(b, x)
def J(self, A, x):
assemble(self.a, tensor=A, reset_sparsity=self.reset_sparsity)
assemble(self.a, tensor=A)
self.bc.apply(A)
self.reset_sparsity = False
# Set parameter values
dt = 0.01 # time step45
T = 1
nu = 0.01
eta = 0.1
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
### Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"
# Define trial and test functions
du = TrialFunction(ME)
(v, q, s) = TestFunctions(ME)
# Define functions
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
dV, dQ, dS = split(du)
V, Q, S = split(u)
V0, Q0, S0 = split(u0)
# Create initial conditions and interpolate
u_init = InitialConditions()
u.interpolate(u_init)
u0.interpolate(u_init)
# Define coefficients
k = Constant(dt)
f = Constant((0, 0))
g = Constant ((0, 9.8))
# Q_(n+theta)
Q_mid = (1.0-theta)*Q0 + theta*Q
# Weak statement of the equations
n = triangle.n
L0 = (1/k)*inner(V-V0, v)*dx + div(V0)*V0*v*dx + nu*inner(grad(V), grad(v))*dx - div(v)*Q_mid*dx- v*dot(g, n)*dx
L1 = div(V)*q*dx
L2 = (1/k)*inner(S, s)*dx-div(V)*S0*s*dx
L = L0 + L1 + L2
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# Create nonlinear problem and Newton solver
problem = NavierStokesEquation(a, L, bc)
newton_solver = NewtonSolver("lu")
newton_solver.parameters["convergence_criterion"] = "incremental"
newton_solver.parameters["relative_tolerance"] = 1e-6
# Output file
file = File("output.pvd", "compressed")
# Step in time
t = 0.0
T = 80*dt
while (t < T):
u0.vector()[:] = u.vector()
newton_solver.solve(problem, u.vector())
t+= dt
file << (u.split()[0], t)
plot(u.split()[0])
interactive()
Could you please help me?
Regards
Saeedeh