This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

The module did not compile for mixed navier-stokes equation

0 votes

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

asked Aug 14, 2014 by saeedeh FEniCS Novice (170 points)
edited Aug 15, 2014 by saeedeh

Please format your code properly, so copy/paste actually produce a working code. The indentation is a mess.

Thanks for your reply. I have edited the code as you suggested. could you please help me in debugging that?

Better, but the indentation is still off. Which version are you using?

Sorry but I copy my codes in code sample exactly. I uses FEniCS 1.0.0. Can I send the module to your email?

There are just a few indents missing. For help, I would suggest upgrading to the latest stable version (1.4.0). 1.0.0 is over three years old.

I think it is OK now. Could you please check it?

...