Hello, thanks for your input! Unfortunately, removing the Dirichlet BC is giving rise to graphs with 0 pressure and velocity values for the entire range of R_in and R_out.

```
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from tqdm import tqdm
plt.rcParams['figure.figsize'] = [14, 8]
plt.rcParams.update({'font.size': 24})
set_log_level(30)
n_elements = 100
mesh = UnitIntervalMesh(n_elements)
# Define a Function Space
U = FunctionSpace(mesh, "Lagrange", 1)
P = FunctionSpace(mesh, "Lagrange", 1)
# Define trial and test functions
u_trial = TrialFunction(U)
p_trial = TrialFunction(P)
u_test = TestFunction(U)
p_test = TestFunction(P)
# The initial condition, u(t=0, x) = sin(pi * x)
u_initial_condition = Expression("0", degree=1)
p_initial_condition = Expression("sin(pi * x[0])", degree=1)
#p_initial_condition = Expression("0", degree=1)
# Discretize the initial condition
u_old = interpolate(u_initial_condition, U)
p_old = interpolate(p_initial_condition, P)
# The time stepping of the implicit Euler discretization (=dt)
EndT = 3 # sec
delta_t = 0.0001
dt = Constant(delta_t)
c = Constant(1)
rho = Constant(1.2)
# The forcing on the rhs of the PDE
heat_source = Constant(0.0)
R_in = 0.5
R_out = 0.5
# Boundary Conditions
#outflow = DirichletBC(P, 0.0, "x[0] > 1 - DOLFIN_EPS")
# inflow = DirichletBC(P, 0.0, "x[0] < DOLFIN_EPS")
p_in = Expression("t < 1.0 ? 0.2*2*3.14*2*sin(2*3.14*2*t) : 0", t=0.0, degree=2)
#inflow = DirichletBC(P, p_in, "x[0] < DOLFIN_EPS")
# bcp = [inflow]
#bcp = [inflow, outflow]
# Define outlet boundary
# Create a MeshFunction for marking boundaries
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
class InletBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0, DOLFIN_EPS) and on_boundary
# Define outlet boundary
class OutletBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0, DOLFIN_EPS) and on_boundary
# Initialize and mark the outlet boundary with a unique ID, e.g., 1
# inlet_boundary = InletBoundary()
# inlet_boundary.mark(boundary_markers, 1)
# outlet_boundary = OutletBoundary()
# outlet_boundary.mark(boundary_markers, 2)
def boundary(x, on_boundary):
return on_boundary
# Define a Measure for the marked outlet boundary
ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)
# Variational formulation
F1 = (inner((p_trial - p_old) / (dt * c * c), p_test) * dx + rho * u_old.dx(0) * p_test * dx
+ (c * (1 + R_in) / (1 - R_in) * p_trial) * p_test * ds(1)
+ (-c * (1 + R_out) / (1 - R_out) * p_trial) * p_test * ds(2))
F2 = inner(((u_trial - u_old) * rho) / dt, u_test) * dx + p_old.dx(0) * u_test * dx
a1 = lhs(F1)
L1 = rhs(F1)
a2 = lhs(F2)
L2 = rhs(F2)
u_solution = Function(U)
p_solution = Function(P)
# Time stepping
n_time_steps = int(np.ceil(EndT / delta_t))
u_plot = np.array(u_old.vector().get_local())
p_plot = np.array(p_old.vector().get_local())
x_coords = mesh.coordinates()[:, 0]
#Boundary points
boundary_indices = np.where((x_coords < DOLFIN_EPS) | (x_coords > 1 - DOLFIN_EPS))[0]
gf_boundary_all = []
time_current = 0.0
for i in tqdm(range(n_time_steps)):
time_current += delta_t
# Update pressure boundary condition
p_in.t = time_current
solve(a1 == L1, p_solution)#, bcp)
solve(a2 == L2, u_solution)
u_old.assign(u_solution)
p_old.assign(p_solution)
u_current = u_solution.vector().get_local()
p_current = p_solution.vector().get_local()
# Compute f and g for the current time step
f = (p_current / (rho(0) * c(0))) + u_current
g = (p_current / (rho(0) * c(0))) - u_current
gf_boundary = np.zeros(len(boundary_indices))
for idx, boundary_index in enumerate(boundary_indices):
if abs(f[boundary_index]) > 1E-12: # Avoid division by zero
gf_boundary[idx] = abs(g[boundary_index] / f[boundary_index])
gf_boundary_all.append(gf_boundary)
u_plot = np.vstack((u_plot, u_current))
p_plot = np.vstack((p_plot, p_current))
gf_boundary_all = np.array(gf_boundary_all)
plt.figure(1)
plt.contourf(u_plot[::100])
plt.colorbar().set_label('$v$ [-]', rotation=270, labelpad=40)
plt.savefig('./Velocity01.png')
plt.figure(2)
plt.contourf(p_plot[::100])
plt.colorbar().set_label('$p$ [-]', rotation=270, labelpad=40)
plt.savefig('./Pressure01.png')
```