Hey!
I want to use the ‘Stokes equations with an iterative solver’ demo as base to calculate with. I was able to change the geometry and boundary conditions to my desire. Now i want to add time to the system an use finite differences just like in the example for the poisson equation with implicit euler.
The altered code is shown below:
from dolfin import *
# Test for PETSc or Tpetra
if not has_linear_algebra_backend("PETSc") and not has_linear_algebra_backend("Tpetra"):
info("DOLFIN has not been configured with Trilinos or PETSc. Exiting.")
exit()
if not has_krylov_solver_preconditioner("amg"):
info("Sorry, this demo is only available when DOLFIN is compiled with AMG "
"preconditioner, Hypre or ML.")
exit()
if has_krylov_solver_method("minres"):
krylov_method = "minres"
elif has_krylov_solver_method("tfqmr"):
krylov_method = "tfqmr"
else:
info("Default linear algebra backend was not compiled with MINRES or TFQMR "
"Krylov subspace method. Terminating.")
exit()
T = 2.0 # final time
num_steps = 500 # number of time steps
dt = T / num_steps # time step size
ufile_pvd = File("Stokes2/Velocitiy/velocity.pvd")
height = 2.0
width = 1.0
length = 1.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length,width,height), 16, 16, 16)
# Build function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
# Boundaries
#Realtion of coordinates:
#x[0] = x = length
#x[1] = y = width
#x[2] = z = height
def Top(x, on_boundary): return x[2] > height - DOLFIN_EPS
#the values for the determination of the boundary must be altered by hand, it is not able to use variables
def Container(x, on_boundary):
return x[0] > length - DOLFIN_EPS or x[0] < DOLFIN_EPS or x[1] > width - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[2] < DOLFIN_EPS
# No-slip boundary condition for velocity on each facet
noslip = Constant((0.0, 0.0, 0.0))
bcContainer = DirichletBC(W.sub(0), noslip, Container)
# Coupling with plate
inflow = Expression(("-sin(x[1]*pi)", "0.0", "0.0"), degree=2)
bcTop = DirichletBC(W.sub(0), inflow, Top)
# Collect boundary conditions
bcs = [bcTop, bcContainer]
#Define variational problem
u_n = Function(W)
phi = TestFunction(W)
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(u,v)*dx + dt*inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(u_n,phi)*dx + dt * inner(f, v)*dx
# Form for use in constructing preconditioner matrix
b = inner(grad(u), grad(v))*dx + p*q*dx
# Create Krylov solver and AMG preconditioner
solver = KrylovSolver(krylov_method, "amg")
U = Function(W)
# Assemble system
A, bb = assemble_system(a, L, bcs)
# Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)
t = 0
for n in range(num_steps):
# Update current time
t += dt
#Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P)
#Solve
solver.solve(U.vector(), bb)
#Get sub-functions
u, p = U.split()
#Save solution in VTK format
ufile_pvd << (u,t)
#Update previous solution
u_n.assign(u)However i get the following error message:
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to initialize vector of degrees of freedom for function.
*** Reason: Cannot re-initialize a non-empty vector. Consider creating a new function.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------I suspect the preconditioner system is the problem and also the assembly of the systems. I couldn’t find any good documentation to work with.
So now the question is how can i successfully add time to my system?
Any help is much appreciated!