TypeError: unsupported operand type(s) for +: 'Sum' and 'Form' Solvi Navier Stokes equation withixed finite elements


82
views
0
3 months ago by
Hi,

I am trying to solve a Navier Stokes PDE problem using FEniCS and mixed elements, but I am facing an error that I do not understand:

TypeError: unsupported operand type(s) for +: 'Sum' and 'Form'

Below is presented my complete code. I hope that someone will be able to help me.

Thank you in advance.

from __future__ import print_function
from fenics import *
from dolfin import *

# Define constants and known data
f = Constant((0.0, 0.0))
mu = Constant(0.005)
rho = Constant(1.0)

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

# Create mesh and define function space
N = 100 # Characteristic number of elements
mesh = RectangleMesh(Point(0,0),Point(1,1),N,N)

# Function space
V = VectorElement('P', mesh.ufl_cell(), 2) # velocity
Q = FiniteElement('P', mesh.ufl_cell(), 1) # pressure
element = MixedElement([V, Q])
W = FunctionSpace(mesh, element)

# Define boundaries 
lid = 'near(x[1], 1.0)'
noslip = 'near(x[1], 0.0) || near(x[0], 0.0) || near(x[0], 1.0)'
pressure = 'near(x[1], 0.0)'

# Boundary conditions
bc_noslip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), noslip)
bc_lid = DirichletBC(W.sub(0), Constant((1.0, 0.0)), lid)
bc_pressure = DirichletBC(W.sub(1), Constant(0.0), pressure)
bc = [bc_noslip, bc_lid, bc_pressure]

# Variational formulation
w = TestFunction(W)
(v, q) = split(w)

# Define the initial conditions and previous time step values
# The previous step's solution
w_n = Function(W)
# Split the mixed function
(u_n, p_n) = split(w_n)
# Create and assign the initial conditions
w_init = Constant((0.0,0.0,0.0)) # null velocity field and pressure 
w_n.interpolate(w_init)
# Create the current solution.
w_ = Function(W)
(u_, p_) = split(w_)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Weak formulation F == 0
Eq1 = dot(u_,v)/dt - dot(u_n,v)/dt \
   + rho*dot(dot(u_, nabla_grad(u_)), v)*dx \
   + 2*mu*inner(epsilon(u_), grad(v))*dx \
   - div(v)*p*dx \
   - dot(f,v)*dx 

Eq2 = q*div(u_)*dx

F = Eq1 + Eq2

# Create the progress bar.
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Provide the Jacobian form.
J = derivative(F,w_)
# Create the Newton solver.
problem = NonlinearVariationalProblem(F, w_, bc, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters.newton_solver
prm.error_on_nonconvergence = False

# Save solution to file in VTK format 
vtkfile1 = File('Navier Stokes Lid Driven Cavity/velocity.pvd')
vtkfile2 = File('Navier Stokes Lid Driven Cavity/pressure.pvd')

# Form and solve the linear system
for n in range(num_steps):

    # Update current time
    t += dt

    # Solve the variational problem
    solver.solve() # for use with the NonlinearVariationalProblem structure
    (u_, p_) = split(w_) # Split solution

    # Save data
    vtkfile1 << (u_, t)
    vtkfile << (p_, t)
    # Update the previous solution
    u_n.assign(u_)

    # Update the progress bar
    progress.update(t/T)​
Community: FEniCS Project

1 Answer


1
3 months ago by
Change the definition of Eq1 to
Eq1 = dot(u_,v)/dt*dx - dot(u_n,v)/dt*dx \
   + rho*dot(dot(u_, nabla_grad(u_)), v)*dx \
   + 2*mu*inner(epsilon(u_), grad(v))*dx \
   - div(v)*p_*dx \
   - dot(f,v)*dx ​
Thank you very much for solving my problem !
written 3 months ago by Matthieu Diaz  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »