Hello, I am trying to simulate a 3D wave equation on a sample box mesh but I have a hard time trying to implement and update boundary conditions.
In my model I have 2 trial functions (u1 is the pressure scalar field and u2 is the velocity vector field); I have 2 equations (first one discretized with implicit euler and the second one with explicit euler), which I summed them up, resulting in this formulation:
\text{(Time-discretized): } \begin{array}{ll} \chi_s\int_\Omega \bar v_1 \cdot \frac{u_1^{n+1}-u_1^n}{\Delta t}\,d\Omega + \mu_0\int_\Omega \bar v_2 \cdot \frac{u_2^{n+1}-u_2^n}{\Delta t}\,d\Omega = \int_\Omega \text{grad}\bar v_1 \cdot \bar u_2^{n+1}\,d\Omega + \int_{\partial \Omega}\bar v_1\cdot \bar u_2^{n+1}\cdot \boldsymbol n\, dS -\int_\Omega \bar v_2\cdot \text{grad}\bar u_1^n\,d\Omega \end{array}
For the functional spaces I chose these ones, and I combined them in order to obtain a mixed space for both pressure and velocity.
P1 = FiniteElement("P-", tetrahedron, 1,0)
N1E = FiniteElement("P-", tetrahedron, 1,1)
Lmix = FunctionSpace(mesh, P1*N1E)
Then I need some functions for the solutions at previous time step (u_n) and for the boundary input (u_b).
# Trial functions
u1, u2 = TrialFunctions(Lmix)
# Test functions
v1, v2 = TestFunctions(Lmix)
# Functions - previous step - boundary inputs
u_n = Function(Lmix) # previous step (with respect to n+1)
u_b = Function(Lmix) # boundary inputs
u_s = Function(Lmix) # solution (to be updated at each time step in the iteration loop)
# Forms
A = v1*u1*dx + dot(v2, u2)*dx - dt*dot(grad(v1), u2)*dx
My problem is in the time stepping loop when I try to define BCs. I tried different syntaxes but none was working properly. This was my working definition for the corresponding 2D problem:
while(t < 0.05):
g = Expression(("0", "sin(pi*t)*near(x[0],0)*(t<1)", "sin(pi*t)*near(x[1],1)*(t<1)"), degree = 2,t = t)
u_b = project(g, Lmix)
K = v1*u_n.split()[0]*dx + dot(v2, u_n.split()[1])*dx - dt*dot(v2, grad(u_n.split()[0]))*dx + \
dt*v1*dot(u_b.split()[1], n)*ds
solve(A==K, u_s, solver_parameters=dict(linear_solver='gmres', preconditioner='ilu'))
# Assign solution to previous step vector
u_n.assign(u_s)
print(t)
# Plot u_1 (pressure)
c = plot(u_s.split()[0])
plt.colorbar(c)
plt.show()
# Update time
t += dt
So for the 3D problem I simply added an additional component (the 3rd velocity) to the Expression definition like that:
g = Expression(("0", "sin(pi*t)*near(x[0],0)*(t<1)", "sin(pi*t)*near(x[1],1)*(t<1)", "0"), degree = 2,t = t)
u_b = project(g, Lmix)
K = v1*u_n.split()[0]*dx + dot(v2, u_n.split()[1])*dx - dt*dot(v2, grad(u_n.split()[0]))*dx + \
dt*v1*dot(u_b.split()[1], n)*ds
But in this case, now it does not return any error but it perform only the first iteration at t=0 and then stops.
I tried to define bcs in different ways, splitting pressure and velocities subspaces like:
Lmix_p = Lmix.sub(0).collapse() # split the mixed functional space to later access the P1 (pressure) subspace
Lmix_v = Lmix.sub(1).collapse() # split the mixed functional space to later access the N1E (velocity) subspace
g_p = Expression('0', degree = 2, t = t)
g_v = Expression(('sin(pi*t)*near(x[0],0)*(t<1)', 'sin(pi*t)*near(x[1],1)*(t<1)', '0'), degree = 2, t = t)
u_b_p = project(g_p, Lmix_p)
u_b_v = project(g_v, Lmix_v)
K = v1*u_n.split()[0]*dx + dot(v2, u_n.split()[1])*dx - dt*dot(v2, grad(u_n.split()[0]))*dx + \
dt*v1*dot(u_b_v, n)*ds
But I still the loop end after one iteration.
Do you have any idea about the proper way to define bcs in order to get the code work in 3D? Maybe I chose worn fucntional spaces?
Any help will be very much appreciated! Thank you very much!
Here is my complete code:
from fenics import *
import numpy as np
# Mesh constants
x0, xL, y0, yL, z0, zL = 0., 0.04, 0., 0.04, 0., 0.05 # meters
z1, z2 = 0.005, 0.045 # interfaces 1 and 2 meters
nx = 10 # divisons x
ny = 10 # divisons y
nz = 25 # divisons z
# Time constants
tinit = 0.
tfinal = 10.
dt = .5
num_steps = (tfinal - tinit)*dt
# Initialize time
t = tinit
# Create mesh and define function space
meshing_start = time.perf_counter()
mesh = BoxMesh(Point(x0, y0, z0), Point(xL, yL, zL), nx, ny, nz)
meshing_time = time.perf_counter()-meshing_start
n = FacetNormal(mesh)
P1 = FiniteElement("P-", tetrahedron, 1,0)
N1E = FiniteElement("P-", tetrahedron, 1,1)
Lmix = FunctionSpace(mesh, P1*N1E)
Lmix_p = Lmix.sub(0).collapse() # split the mixed functional space to later access the P1 (pressure) subspace
Lmix_v = Lmix.sub(1).collapse() # split the mixed functional space to later access the N1E (velocity) subspace
# Define boundaries
def boundary(x):
return near(x[0],0) or near(x[0],xL) or near(x[1],0) or near(x[1],yL)or near(x[2],0) or near(x[2],zL)
# Trial functions
u1, u2 = TrialFunctions(Lmix)
# Test functions
v1, v2 = TestFunctions(Lmix)
# Functions - previous step - boundary inputs
u_n = Function(Lmix) # previous step (with respect to n+1)
u_b = Function(Lmix) # boundary inputs
u_s = Function(Lmix) # solution (to be updated at each time step in the iteration loop)
# Forms
A = v1*u1*dx + dot(v2, u2)*dx - dt*dot(grad(v1), u2)*dx
while(t < 0.5):
g = Expression(("0", "sin(pi*t)*near(x[0],0)*(t<1)", "sin(pi*t)*near(x[1],1)*(t<1)", "0"), degree = 2,t = t)
u_b = project(g, Lmix)
K = v1*u_n.split()[0]*dx + dot(v2, u_n.split()[1])*dx - dt*dot(v2, grad(u_n.split()[0]))*dx + \
dt*v1*dot(u_b.split()[1], n)*ds
solve(A==K, u_s, solver_parameters=dict(linear_solver='gmres', preconditioner='ilu'))
# Assign solution to previous step vector
u_n.assign(u_s)
print(t)
# Plot u_1 (pressure)
c = plot(u_s.split()[0])
plt.colorbar(c)
plt.show()
# Update time
t += dt
Niccolò