Good evening,
I am trying to solve the unsteady Navier Stokes equation in a complex geometry, and for now everything worked fine using Oasisx.
However, I want to try to implement a time-dependent Dirichlet boundary condition for the pressure, but I cannot implement it.
As my code is very long, I implemented a time-dependent pressure condition starting from this code and also referring to the definition of oasisx.PressureBC here.
The MWE is:
from mpi4py import MPI
import dolfinx
import numpy as np
import oasisx
from oasisx import DirichletBC, LocatorMethod, FractionalStep_AB_CN
from typing import List
def inlet(x):
return np.isclose(x[0], 0)
def wall(x):
return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))
def outlet(x):
return np.isclose(x[0], 10)
domain = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([10, 1])], [
50, 5], cell_type=dolfinx.mesh.CellType.triangle)
bcs_p: List[oasisx.PressureBC] = []
bc_inlet_x = DirichletBC(dolfinx.fem.Constant(domain, 1.),
method=LocatorMethod.GEOMETRICAL, marker=inlet)
bc_inlet_y = DirichletBC(dolfinx.fem.Constant(domain, 0.),
method=LocatorMethod.GEOMETRICAL, marker=inlet)
bc_wall = DirichletBC(dolfinx.fem.Constant(domain, 0.), method=LocatorMethod.GEOMETRICAL,
marker=wall)
outlet_facets = dolfinx.mesh.locate_entities_boundary(domain, domain.topology.dim - 1, outlet)
ft = dolfinx.mesh.meshtags(domain, domain.topology.dim-1, outlet_facets,
np.full_like(outlet_facets, 1, dtype=np.int32))
class Pressure_Outlet():
def __init__(self, t:float):
self.t = t
def eval(self,x: np.typing.NDArray[np.float64]):
return 1*self.t
t = 0.
p_val = Pressure_Outlet(t)
bcs_p: List[oasisx.PressureBC] = [oasisx.PressureBC(p_val.eval, (ft, 1))]
bcs_u = [[bc_inlet_x, bc_wall], [bc_wall, bc_inlet_y]]
# fractional step solver
solver = FractionalStep_AB_CN(
mesh=domain,
u_element=("Lagrange", 2),
p_element=("Lagrange", 1),
bcs_u=bcs_u,
bcs_p=bcs_p,
solver_options={"tentative": {"ksp_type": "preonly", "pc_type": "lu"}, "pressure": {
"ksp_type": "preonly", "pc_type": "lu"}, "scalar": {"ksp_type": "preonly", "pc_type": "lu"}},
body_force=None
)
# Time-stepping
T_start, T_end, dt = 0.0, 1, 0.001
num_steps = int((T_end - T_start) // dt)
with dolfinx.io.VTXWriter(domain.comm, "poisseuille.bp", [solver.u], engine="BP4") as writer:
for step in range(num_steps):
t += dt
p_val.t = t
solver.solve(dt, nu=0.01, max_iter=10)
writer.write(dt*step)
And the error is:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
File ~/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/oasisx/bcs.py:215, in PressureBC.create_bcs(self, V, Q)
214 try:
--> 215 rhs = [self._value * n_i * v.dx(i) * ds for i, n_i in enumerate(n)]
216 except TypeError:
217 # If input is lambda function interpolate into local function
TypeError: unsupported operand type(s) for *: 'method' and 'Indexed'
During handling of the above exception, another exception occurred:
TypeError Traceback (most recent call last)
File ~/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py:397, in Function.interpolate(self, u, cells, nmm_interpolation_data)
395 try:
396 # u is a Function or Expression (or pointer to one)
--> 397 _interpolate(u, cells)
398 except TypeError:
399 # u is callable
File ~/anaconda3/envs/fenicsx-env/lib/python3.12/functools.py:909, in singledispatch.<locals>.wrapper(*args, **kw)
906 raise TypeError(f'{funcname} requires at least '
907 '1 positional argument')
--> 909 return dispatch(args[0].__class__)(*args, **kw)
File ~/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py:373, in Function.interpolate.<locals>._interpolate(u, cells)
372 """Interpolate a cpp.fem.Function"""
--> 373 self._cpp_object.interpolate(u, cells, nmm_interpolation_data)
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
1. (self: dolfinx.cpp.fem.Function_float64, f: numpy.ndarray[numpy.float64], cells: numpy.ndarray[numpy.int32]) -> None
2. (self: dolfinx.cpp.fem.Function_float64, u: dolfinx.cpp.fem.Function_float64, cells: numpy.ndarray[numpy.int32], nmm_interpolation_data: Tuple[List[int], List[int], List[float], List[int]]) -> None
3. (self: dolfinx.cpp.fem.Function_float64, expr: dolfinx::fem::Expression<double, double>, cells: numpy.ndarray[numpy.int32]) -> None
Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7f3b3720d0b0>, <bound method Pressure_Outlet.eval of <__main__.Pressure_Outlet object at 0x7f3b372b4110>>, array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155,
156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168,
169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181,
182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194,
195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,
208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220,
221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233,
234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246,
247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259,
260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272,
273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285,
286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298,
299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311,
312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324,
325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337,
338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350,
351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363,
364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376,
377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389,
390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402,
403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415,
416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428,
429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441,
442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454,
455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467,
468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480,
481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493,
494, 495, 496, 497, 498, 499], dtype=int32), ((), (), (), ())
During handling of the above exception, another exception occurred:
IndexError Traceback (most recent call last)
Cell In[3], line 49
44 bcs_u = [[bc_inlet_x, bc_wall], [bc_wall, bc_inlet_y]]
48 # fractional step solver
---> 49 solver = FractionalStep_AB_CN(
50 mesh=domain,
51 u_element=("Lagrange", 2),
52 p_element=("Lagrange", 1),
53 bcs_u=bcs_u,
54 bcs_p=bcs_p,
55 solver_options={"tentative": {"ksp_type": "preonly", "pc_type": "lu"}, "pressure": {
56 "ksp_type": "preonly", "pc_type": "lu"}, "scalar": {"ksp_type": "preonly", "pc_type": "lu"}},
57 body_force=None
58 )
61 # Time-stepping
62 T_start, T_end, dt = 0.0, 1, 0.001
File ~/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/oasisx/fracstep.py:198, in FractionalStep_AB_CN.__init__(self, mesh, u_element, p_element, bcs_u, bcs_p, solver_options, jit_options, body_force, options)
196 self._bcs_p = bcs_p
197 for bcp in self._bcs_p:
--> 198 bcp.create_bcs(self._Vi[0][0], self._Q)
199 for i in range(self._mesh.geometry.dim):
200 forms[i].append(bcp.rhs(i))
File ~/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/oasisx/bcs.py:219, in PressureBC.create_bcs(self, V, Q)
216 except TypeError:
217 # If input is lambda function interpolate into local function
218 self._u = _fem.Function(Q)
--> 219 self._u.interpolate(self._value) # type: ignore
220 rhs = [self._u * n_i * v.dx(i) * ds for i, n_i in enumerate(n)]
222 # Create rhs contribution from natural boundary condition
File ~/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py:402, in Function.interpolate(self, u, cells, nmm_interpolation_data)
400 assert callable(u)
401 x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh.geometry, cells)
--> 402 self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)
IndexError: invalid axis: 0 (ndim = 0)
I do not understand what am I doing wrong, as I implemented other time-dependent Dirichlet boundary condition for the velocity with a similar sintax without getting any error.
Then, as I am asking about time dependent pressure boundary condition, I know that using this kind of boundary condition with IPCS scheme is not well supported, and I tried to follow the suggestions given in this reply modifying the scheme as suggested, however even if I got overall convergence I have unphysical solution on the boundary where the dirichlet pressure is imposed.
Any suggestion to which scheme could better suit a time-dependent pressure profile to be imposed on a boundary?
Thank you for your time.