Hi, quick question.
I want to calculate a coupled system and i need the derivative of the solution of a previous calculated variational problem as the boundary condition for the next equation. Is there a fast way to save the solution in a 2D array, so i can build the derivative and use as a boundary condition?
Thanks in advance!
Please add a minimal working code example.
You need to clarify what the derivative of the solution is with respect to (is it a spatial derivative, or something else).
Hi Dokken,
I want to couple the biharmonic equation with the stokes equation. The displacement and the following derivative in time of the plate which is modelled by the biharmonic equation is supposed to be the boundary condition of the fluid system which is modelled by the stokes equation.
Following code produces quite good results for a evolution in time:
mesh = UnitSquareMesh(16,16)
X = FiniteElement('CG', mesh.ufl_cell(), 2)
Y = FiniteElement('CG', mesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(mesh, Z)
u = TrialFunctions(S)
phi = TestFunctions(S)
a = inner(grad(u[1]), grad(phi[0]))*dx + inner(grad(u[0]), grad(phi[1]))*dx
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(S, Constant((0,0)), boundary)
L_dummy = Constant(0)*phi[0]*dx
m = u[0]*phi[0]*dx + u[1]*phi[1]*dx
A, _ = assemble_system(a,L_dummy,bc)
B = assemble(m)
eigensolver = SLEPcEigenSolver(as_backend_type(A), as_backend_type(B))
eigensolver.parameters['spectrum'] = 'target real'
eigensolver.parameters['tolerance'] = 1e-6
eigensolver.solve()
r, c, rx, cx = eigensolver.get_eigenpair(0)
print("Eigenvalue: ", r)
ef = Function(S)
ef.vector()[:] = rx
#plot(ef[0])
#
# %% Computation of time dependent problem
#
T = 1.5 # final time
num_steps = 5000 # number of time steps
dt = T / num_steps # time step size
t = 0
u_m = Function(S)
u_n = Function(S)
u = TrialFunctions(S)
phi = TestFunctions(S)
# Create VTK file for saving solution
vtkfile = File('solution/ResonatingPlate.pvd')
#define F and assemble unknown and known variables
a_evo = -dt*dt*(inner(grad(u[0]),grad(phi[0])) + inner(grad(u[1]),grad(phi[1])) + u[1]*phi[0])*dx + (u[0]*phi[1])*dx
f = Constant(0.0)
L = (2*u_n[0]*phi[1] - u_m[0]*phi[1])*dx + dt*dt*(inner(f*ef[0],phi[0]))*dx
l=Function(S)
for n in range(num_steps):
# Compute solution
solve(a_evo == L, l, bc)
# Save to file
vtkfile << (l, t)
#vtkfile2 << (w_temp,t)
u_m.assign(u_n)
# Update previous solution
u_n.assign(l)
# Update current time
t += dt
f = Constant(2*r*sin(r*t))
L = (2*u_n[0]*phi[1] - u_m[0]*phi[1])*dx + dt*dt*(inner(f*ef[0],phi[0]))*dx
# Update time of the expression
f.t = tNow i need to find a way of using the results in the following codesection for the stokes system:
# 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()
height = 1.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)
def Top(x, on_boundary): return x[2] > height - DOLFIN_EPS
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)Here it gets tricky, i basically want to use use the results as a cover wihich stimulates the fluid, the resulting pressure that works back at the plate is an additional task. Moreover i am right now not sure how to handle the preconditioner matrices.
# Coupling with plate
#inflow = Expression(('0.0, 0.0, displacement*x[2]'), degree = 2)
bcTop = DirichletBC(W.sub(0), (0.0, 0.0, displacement), Top)
# Collect boundary conditions
bcs = [bcTop, bcContainer]
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx
# Form for use in constructing preconditioner matrix
b = inner(grad(u), grad(v))*dx + p*q*dx
# Assemble system
A, bb = assemble_system(a, L, bcs)
# Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)
# Create Krylov solver and AMG preconditioner
solver = KrylovSolver(krylov_method, "amg")
# Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P)
# Solve
U = Function(W)
solver.solve(U.vector(), bb)
# Get sub-functions
u, p = U.split()
# Save solution in VTK format
ufile_pvd = File("CoupledSystem/Velocitiy/velocity.pvd")
ufile_pvd << u
pfile_pvd = File("CoupledSystem/Pressure/pressure.pvd")
pfile_pvd << p Hi, i found a way to build the derivative, so that point is clear now. The solution was a simple differential quotient.
Now i need to use the derivative as a boundary condition in every time step. After digging up some older posts i found an answer of you, how to save the solution of a function in an array. My problem is now that the boundary condition expects a 3-dimensional input vor the values.
Is there any other way to use calculated values as a boundary condition? I thought about just using the maximum value of the resulting velocity profile (derivative), but this seems too simple and not really expedient.
Thanks in advance
Please simplify your code such that it is clear where the issue is. You do not need to solve any actual physical problem to illustrate your problem. Define the variables you would have had in the original problem, and what their relation is. Having the mathematical formulation of what you would like to achieve would also help.
Okay, sorry that i wasn’t that clear, what i want.
My main problem is that i calculate a function over an 2-dimensional area, which was calculated in the first code snippet.
The solution itself is a function of the following function space:
X = FiniteElement('CG', mesh.ufl_cell(), 2)
Y = FiniteElement('CG', mesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(mesh, Z)
l=Function(S)This solution needs to be derived, that point is no problem anymore.
Now the real problem is to use the calculated function l as input for the following boundary condition:
def Top(x, on_boundary): return x[2] > height - DOLFIN_EPS
bcTop = DirichletBC(W.sub(0), (0.0, 0.0, PlateSpeed), Top)Here arises the problem because i am not sure how get the values of the calculated function and plug them into the Dirichlet bc. Though it is correct that i only need the values in z-axis because the plate only moves in z-direction i want to insert all the values of the function l over the area that is the top cover of the volume which is modelled by the stokes equation.
I added a sketch with a formulation and a visualization from my professor.
If plate-speed is (phi-phi_n)/dt where phi is the solution at the current time step and phi_n the solution from the previous time step, you could project this into an appropriate function space to apply the boundary condition.
W0=W.sub(0).collapse()
bc_function =project(as_vector((0,0, (phi-phi_n)/dt)),W0)
and then use this in your boundary condition.
Hi, i tried to use your code, however i get an error message indicating, that there are some parts not implemented yet:
Exception: codim != 0 or 1 - Not (yet) implementedThe version i am using is: 2019.2.0.dev0
You ned to produce a minimal code example that reproduces the error (that I can run on my own computer) for me to be able to help you.
My first try to test your hint was to use a placeholder solution, but from the same function space i will work with later.
This is the following snippet:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(16,16)
X = FiniteElement('CG', mesh.ufl_cell(), 2)
Y = FiniteElement('CG', mesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(mesh, Z)
u = TrialFunctions(S)
phi = TestFunctions(S)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(S, Constant((0,0)), boundary)
class Source(UserExpression):
def eval(self, values, x):
values[0] = 4.0*pi**4*sin(pi*x[0])*sin(pi*x[1])
a = inner(grad(u[1]), grad(phi[0]))*dx + inner(grad(u[0]), grad(phi[1]))*dx
f = Source(degree=2)
# Define linear form
L = inner(f,phi[1])*dx
# Solve variational problem
u_Platte = Function(S)
solve(a == L, u_Platte, bc)
# Plot solution
plot(u_Platte[0])
dfile_pvd = File("Stokes2/Plate/Displacement.pvd")
dfile_pvd << u_PlatteThis works perfectly fine.
The next snippet cuases an error in the last line:
#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()
height = 1.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)
def Top(x, on_boundary): return x[2] > height - DOLFIN_EPS
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 side
noslip = Constant((0.0, 0.0, 0.0))
bcContainer = DirichletBC(W.sub(0), noslip, Container)
# Coupling with plate
W0 = W.sub(0).collapse()
placeholder = u_Platte[0]
bc_function = project(as_vector((0,0,placeholder)),W0)Please state what your error message.Please also remove all code not required to produce the error, including variables and if-tests.
The error message is:
Exception: codim != 0 or 1 - Not (yet) implementedThe code, that produces the error is:
W0 = W.sub(0).collapse()
bc_function = project(as_vector((0,0,u_Platte[0])),W0)Code that is lonked to last snippet is:
mesh = UnitSquareMesh(16,16)
X = FiniteElement('CG', mesh.ufl_cell(), 2)
Y = FiniteElement('CG', mesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(mesh, Z)
u_Platte = Function(S)and
height = 1.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)I get no errors running the following:
from dolfin import *
mesh = UnitSquareMesh(10,10)
X = FiniteElement('CG', mesh.ufl_cell(), 2)
Y = FiniteElement('CG', mesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(mesh, Z)
u_Platte = Function(S)
height = 1.0
width = 1.0
length = 1.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length,width,height), 10,10,10)
# Build function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
W0 = W.sub(0).collapse()
bc_function = project(as_vector((0,0,u_Platte[0])),W0, solver_type="mumps")
I tried the addition to the last line, but i still get the same error
What versjon of dolfin are you using?
The version i am using is: 2019.2.0.dev0
I can reproduce the issue using the development version (2019.2.0.dev0). However, could you clarify a few things for me.
Why are you projecting a quantity from a 2D mesh to a 3D mesh? Shouldn’t these quantities live on the same mesh?
What version did you use to get no errors?
Yes, of course: well the biharmonic equation couldnt be solved with H2 elements and using the method to solve it like in the tutorial wasn’t helpfull either, because we also need eigenvalues and eigenfunctions and those were far from correct. So we split the system and worked on the biharmonic equation, now we want to couple a resonating plate with a fluid. So the reason for the different function spaces stems from the two systems being viewed indivudial.
We first wanted to have a plate, that can be put i a resonating state, which works very good. Now we want to see wether an fluid can damp that the resonance
For me the code ran with 2019.1.0.
I think you should consider using MeshView’s, as introduced in
With examples:
This introduced 3D-2D coupling of PDEs
Thanks for the input. I got the an eigenvalue and eigenfunction calculation running as well as calculation for some displacement, but i still cant solve my issue. I build a submesh the following way:
height = 1.0
width = 1.0
length = 1.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length,width,height), 12, 12, 12)
plot(mesh)
#%%
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for f in facets(mesh):
marker[f] = height - DOLFIN_EPS < f.midpoint().z() < height + DOLFIN_EPS
## Create submesh ##
submesh = MeshView.create(marker, 1)The boxmesh will be used for the fluid as i mentioned.
I got an function from the submesh:
X = FiniteElement('CG', submesh.ufl_cell(), 2)
Y = FiniteElement('CG', submesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(submesh, Z)
u_Platte = Function(S)intuitively i would insert the function the following way:
def Top(x, on_boundary): return x[2] > height - DOLFIN_EPS
bcTop = DirichletBC(W.sub(0), u_Platte[0], Top)But i get the following error:
UFLException: Shapes do not match: <Argument id=140260893502624> and <Indexed id=140261033632928>.How do i math the shapes? The function maybe from different function space, but from the same mesh.