Hey guys,
below you find my attempt to implement the navier stokes eq. step by step. now i am just trying to solve two of t hose and work my way up:
from __future__ import print_function
import fenics
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from mshr import *
from fenics import XDMFFile
from fenics import derivative
# Mesh und Umgebung erstellen
l1 = 3
l2 = 1
domain = Rectangle(Point(0, 0), Point(l1, l2))
mesh = generate_mesh(domain, 128)
# Elemente, Funktionen und Funktionsraum erstellen
P2 = VectorElement("P", triangle, 2)
P1 = FiniteElement("P", triangle, 1)
element = MixedElement([P2, P1])
V = FunctionSpace(mesh, element)
solution = Function(V)
u, rho = solution.split()
solution_n = Function(V)
u_n, rho_n = solution_n.split()
w_1, w_2 = TestFunctions(V)
# Zeitdiskretisierung
T = 5.0 # final time
f = 1/T
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
# erste vereinfachte Gleichung
F1 = inner((u-u_n), w_1) * dx \
+ dt * inner(nabla_grad(u), nabla_grad(w_1)) * dx
F2 = inner((rho-rho_n), w_2) * dx \
+ dt * inner(div(u), w_2) * dx
F = F1 + F2
a = lhs(F)
L = rhs(F)
# Anfangsbedingungen
u_init = Constant((0.0, 0.0))
u_1_init = Constant(0.0)
u_2_init = Constant(0.0)
rho_init = Constant(998.0)
solution_init = Expression(('u_1_init', 'u_2_init', 'rho_init'), u_1_init = u_1_init, u_2_init = u_2_init, rho_init = rho_init, degree=2)
# Randbedingungen + Schallwelle
t = 0
amplitudeu = Constant((1.5, 0.0))
Schallwelle = Expression(('u_init[0] + amplitudeu[0] * sin(2 * pi * f * t)', 'u_init[1] + amplitudeu[1] * sin(2 * pi * f * t)'), u_init=u_init, amplitudeu=amplitudeu, f=f, t=t, degree=2, domain=mesh)
Wand = Constant((0.0, 0.0))
tol = 1E-14
def boundary_Dl(x, on_boundary):
if on_boundary:
if near(x[0], 0, tol):
return True
else:
return False
else:
return False
def boundary_Dr(x, on_boundary):
if on_boundary:
if near(x[0], l1, tol):
return True
else:
return False
else:
return False
bcl = DirichletBC(V.sub(0), Schallwelle, boundary_Dl)
bcr = DirichletBC(V.sub(0), Wand, boundary_Dr)
bc = [bcl, bcr]
# pvd files erstellen´
xdmf_file_u = XDMFFile('RAWHEATNAVIERtest/u_1.xdmf')
xdmf_file_rho = XDMFFile('RAWHEATNAVIERtest/rho_1.xdmf')
solutionX = Function(V)
solution = Function(V) **LINE 100**
t = 0
for n in range(num_steps):
# Solve for u
solve(a == L, solution, bc)
# Assign values to the storage function
solutionX.assign(solution)
# Write to XDMF file
xdmf_file_u.write(u_solution, t)
# Update for the next time step
t += float(dt)
Schallwelle.t = t
bc = [bcl, bcr]
u_n.assign(u)
This leads to the following error:
Form has no parts with arity 2.
Traceback (most recent call last):
File “RAWHEATNAVIER.py”, line 115, in
solve(a == L, solution, bc)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 220, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 242, in _solve_varproblem
form_compiler_parameters=form_compiler_parameters)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/problem.py”, line 59, in init
cpp.fem.LinearVariationalProblem.init(self, a, L, u._cpp_object, bcs)
RuntimeError:
*** Error: Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason: Expecting the left-hand side to be a bilinear form (not rank 0).
*** Where: This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0
I read that the problem might be the mixed element space but I dont really know what to do differently. Thank you very much for reading and helping!