Skip to content

Commit 2f59aca

Browse files
simulkadesimulkade
authored
Merge pull request #71 from mhvwerts/main
Fix left/right periodic boundary conditions in Grid3D. Block use of "radial periodic BCs". Fix top/bottom BCs in 2D.
2 parents 3427e1d + 91eabb0 commit 2f59aca

4 files changed

Lines changed: 355 additions & 152 deletions

File tree

‎src/pyfvtool/advection.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1025,8 +1025,8 @@ def convectionUpwindTerm3D(u: FaceVariable, *args):
10251025
AN[:, -1, :] = AN[:, -1, :]/2.0
10261026
APy[:, -1, :] = APy[:, -1, :]+vn_min[:, -1, :]/(2.0*DYp[:, -1, :])
10271027
# Back boundary:
1028-
APz[:, :, 0] = APz[:, :, 1]-wb_max[:, :, 1]/(2.0*DZp[:, :, 0])
1029-
AB[:, :, 0] = AB[:, :, 1]/2.0
1028+
APz[:, :, 0] = APz[:, :, 0]-wb_max[:, :, 0]/(2.0*DZp[:, :, 0])
1029+
AB[:, :, 0] = AB[:, :, 0]/2.0
10301030
# Front boundary:
10311031
AF[:, :, -1] = AF[:, :, -1]/2.0
10321032
APz[:, :, -1] = APz[:, :, -1]+wf_min[:, :, -1]/(2.0*DZp[:, :, -1])

‎src/pyfvtool/boundary.py‎

Lines changed: 162 additions & 133 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66

77
from .mesh import MeshStructure
88
from .mesh import Grid1D, Grid2D, Grid3D
9+
from .mesh import SphericalGrid1D, CylindricalGrid1D
910
from .mesh import CylindricalGrid2D
1011
from .mesh import PolarGrid2D, CylindricalGrid3D, SphericalGrid3D
1112
from .utilities import int_range
@@ -897,6 +898,9 @@ def boundaryConditionsTerm1D(BC: BoundaryConditions1D):
897898
s[q] = -(BC.left.b.item()/2 - BC.left.a.item()/dx_1)
898899
BCRHS[G[i]] = -BC.left.c.item()
899900
elif BC.right.periodic or BC.left.periodic: # periodic boundary condition
901+
if (type(BC.domain) is SphericalGrid1D)\
902+
or (type(BC.domain) is CylindricalGrid1D):
903+
raise ValueError("Radial periodic boundary conditions are not physically meaningful.")
900904
# Right boundary
901905
i = Nx+1
902906
q = q+1
@@ -1024,11 +1028,11 @@ def boundaryConditionsTerm2D(BC: BoundaryConditions2D):
10241028
s[q] = 1
10251029
q = q[-1]+i
10261030
ii[q] = G[i,j]
1027-
jj[q] = G[i,Ny+1]
1031+
jj[q] = G[i,Ny]
10281032
s[q] = -1
10291033
q = q[-1]+i
10301034
ii[q] = G[i,j]
1031-
jj[q] = G[i,Ny+2]
1035+
jj[q] = G[i,Ny+1]
10321036
s[q] = -1
10331037
BCRHS[G[i,j]] = 0
10341038

@@ -1058,6 +1062,8 @@ def boundaryConditionsTerm2D(BC: BoundaryConditions2D):
10581062
s[q] = -(BC.left.b/2 - BC.left.a/dx_1)
10591063
BCRHS[G[i,j]] = -BC.left.c
10601064
elif BC.right.periodic or BC.left.periodic: # periodic boundary condition
1065+
if (type(BC.domain) is CylindricalGrid2D):
1066+
raise ValueError("Radial periodic boundary conditions are not physically meaningful.")
10611067
# Right boundary
10621068
i = Nx+1
10631069
j = int_range(1, Ny)
@@ -1262,7 +1268,7 @@ def boundaryConditionsTerm3D(BC: BoundaryConditions3D):
12621268
ii[q] = G[i,j,k].ravel()
12631269
jj[q] = G[0,j,k].ravel()
12641270
s[q] = dx_end/dx_1
1265-
q = q[-1]+int_range[1,Ny*Nz]
1271+
q = q[-1]+int_range(1,Ny*Nz)
12661272
ii[q] = G[i,j,k].ravel()
12671273
jj[q] = G[1,j,k].ravel()
12681274
s[q] = -dx_end/dx_1
@@ -1454,11 +1460,11 @@ def boundaryConditionsTermPolar2D(BC: BoundaryConditions2D):
14541460
s[q] = 1
14551461
q = q[-1]+i
14561462
ii[q] = G[i,j]
1457-
jj[q] = G[i,Ny+1]
1463+
jj[q] = G[i,Ny]
14581464
s[q] = -1
14591465
q = q[-1]+i
14601466
ii[q] = G[i,j]
1461-
jj[q] = G[i,Ny+2]
1467+
jj[q] = G[i,Ny+1]
14621468
s[q] = -1
14631469
BCRHS[G[i,j]] = 0
14641470

@@ -1488,45 +1494,53 @@ def boundaryConditionsTermPolar2D(BC: BoundaryConditions2D):
14881494
s[q] = -(BC.left.b/2 - BC.left.a/dx_1)
14891495
BCRHS[G[i,j]] = -BC.left.c
14901496
elif BC.right.periodic or BC.left.periodic: # periodic boundary condition
1491-
# Right boundary
1492-
i = Nx+1
1493-
j = int_range(1, Ny)
1494-
q = q[-1]+j
1495-
ii[q] = G[i,j]
1496-
jj[q] = G[i,j]
1497-
s[q] = 1
1498-
q = q[-1]+j
1499-
ii[q] = G[i,j]
1500-
jj[q] = G[i-1,j]
1501-
s[q] = -1
1502-
q = q[-1]+j
1503-
ii[q] = G[i,j]
1504-
jj[q] = G[0,j]
1505-
s[q] = dx_end/dx_1
1506-
q = q[-1]+j
1507-
ii[q] = G[i,j]
1508-
jj[q] = G[1,j]
1509-
s[q] = -dx_end/dx_1
1510-
BCRHS[G[i,j]] = 0
1511-
# Left boundary
1512-
i = 0
1513-
q = q[-1]+j
1514-
ii[q] = G[i,j]
1515-
jj[q] = G[i,j]
1516-
s[q] = 1.0
1517-
q = q[-1]+j
1518-
ii[q] = G[i,j]
1519-
jj[q] = G[i+1,j]
1520-
s[q] = 1.0
1521-
q = q[-1]+j
1522-
ii[q] = G[i,j]
1523-
jj[q] = G[Nx,j]
1524-
s[q] = -1.0
1525-
q = q[-1]+j
1526-
ii[q] = G[i,j]
1527-
jj[q] = G[Nx+1,j]
1528-
s[q] = -1.0
1529-
BCRHS[G[i,j]] = 0.0
1497+
raise ValueError("Radial periodic boundary conditions are not physically meaningful.")
1498+
#
1499+
# Keep the following code for future reference, once a physically relevant
1500+
# case has been identified for radial periodic BCs...
1501+
#
1502+
# # Right boundary
1503+
# i = Nx+1
1504+
# j = int_range(1, Ny)
1505+
# q = q[-1]+j
1506+
# ii[q] = G[i,j]
1507+
# jj[q] = G[i,j]
1508+
# s[q] = 1
1509+
# q = q[-1]+j
1510+
# ii[q] = G[i,j]
1511+
# jj[q] = G[i-1,j]
1512+
# s[q] = -1
1513+
# q = q[-1]+j
1514+
# ii[q] = G[i,j]
1515+
# jj[q] = G[0,j]
1516+
# s[q] = dx_end/dx_1
1517+
# q = q[-1]+j
1518+
# ii[q] = G[i,j]
1519+
# jj[q] = G[1,j]
1520+
# s[q] = -dx_end/dx_1
1521+
# BCRHS[G[i,j]] = 0
1522+
# # Left boundary
1523+
# i = 0
1524+
# q = q[-1]+j
1525+
# ii[q] = G[i,j]
1526+
# jj[q] = G[i,j]
1527+
# s[q] = 1.0
1528+
# q = q[-1]+j
1529+
# ii[q] = G[i,j]
1530+
# jj[q] = G[i+1,j]
1531+
# s[q] = 1.0
1532+
# q = q[-1]+j
1533+
# ii[q] = G[i,j]
1534+
# jj[q] = G[Nx,j]
1535+
# s[q] = -1.0
1536+
# q = q[-1]+j
1537+
# ii[q] = G[i,j]
1538+
# jj[q] = G[Nx+1,j]
1539+
# s[q] = -1.0
1540+
# BCRHS[G[i,j]] = 0.0
1541+
#
1542+
#
1543+
#
15301544
# Build the sparse matrix of the boundary conditions
15311545
q = q[-1] + 1
15321546
BCMatrix = csr_array((s[0:q], (ii[0:q], jj[0:q])),
@@ -1677,49 +1691,57 @@ def boundaryConditionsTermCylindrical3D(BC: BoundaryConditions3D):
16771691
s[q] = -(BC.left.b/2 - BC.left.a/dx_1).ravel()
16781692
BCRHS[G[i,j,k].ravel()] = -(BC.left.c).ravel()
16791693
elif BC.right.periodic or BC.left.periodic: # periodic
1680-
# Right boundary
1681-
i=Nx+1
1682-
j=j_ind
1683-
k=k_ind
1684-
q = q[-1]+int_range(1,Ny*Nz)
1685-
ii[q] = G[i,j,k].ravel()
1686-
jj[q] = G[i,j,k].ravel()
1687-
s[q] = 1.0
1688-
q = q[-1]+int_range(1,Ny*Nz)
1689-
ii[q] = G[i,j,k].ravel()
1690-
jj[q] = G[i-1,j,k].ravel()
1691-
s[q] = -1.0
1692-
q = q[-1]+int_range(1,Ny*Nz)
1693-
ii[q] = G[i,j,k].ravel()
1694-
jj[q] = G[0,j,k].ravel()
1695-
s[q] = dx_end/dx_1
1696-
q = q[-1]+int_range[1,Ny*Nz]
1697-
ii[q] = G[i,j,k].ravel()
1698-
jj[q] = G[1,j,k].ravel()
1699-
s[q] = -dx_end/dx_1
1700-
BCRHS[G[i,j,k].ravel()] = 0.0
1701-
1702-
# Left boundary
1703-
i = 0
1704-
j=j_ind
1705-
k=k_ind
1706-
q = q[-1]+int_range(1,Ny*Nz)
1707-
ii[q] = G[i,j,k].ravel()
1708-
jj[q] = G[i,j,k].ravel()
1709-
s[q] = 1.0
1710-
q = q[-1]+int_range(1,Ny*Nz)
1711-
ii[q] = G[i,j,k].ravel()
1712-
jj[q] = G[i+1,j,k].ravel()
1713-
s[q] = 1.0
1714-
q = q[-1]+int_range(1,Ny*Nz)
1715-
ii[q] = G[i,j,k].ravel()
1716-
jj[q] = G[Nx,j,k].ravel()
1717-
s[q] = -1.0
1718-
q = q[-1]+int_range(1,Ny*Nz)
1719-
ii[q] = G[i,j,k].ravel()
1720-
jj[q] = G[Nx+1,j,k].ravel()
1721-
s[q] = -1.0
1722-
BCRHS[G[i,j,k].ravel()] = 0.0
1694+
raise ValueError("Radial periodic boundary conditions are not physically meaningful.")
1695+
#
1696+
# Keep the following code for future reference, once a physically relevant
1697+
# case has been identified for radial periodic BCs...
1698+
#
1699+
# # Right boundary
1700+
# i=Nx+1
1701+
# j=j_ind
1702+
# k=k_ind
1703+
# q = q[-1]+int_range(1,Ny*Nz)
1704+
# ii[q] = G[i,j,k].ravel()
1705+
# jj[q] = G[i,j,k].ravel()
1706+
# s[q] = 1.0
1707+
# q = q[-1]+int_range(1,Ny*Nz)
1708+
# ii[q] = G[i,j,k].ravel()
1709+
# jj[q] = G[i-1,j,k].ravel()
1710+
# s[q] = -1.0
1711+
# q = q[-1]+int_range(1,Ny*Nz)
1712+
# ii[q] = G[i,j,k].ravel()
1713+
# jj[q] = G[0,j,k].ravel()
1714+
# s[q] = dx_end/dx_1
1715+
# q = q[-1]+int_range(1,Ny*Nz)
1716+
# ii[q] = G[i,j,k].ravel()
1717+
# jj[q] = G[1,j,k].ravel()
1718+
# s[q] = -dx_end/dx_1
1719+
# BCRHS[G[i,j,k].ravel()] = 0.0
1720+
1721+
# # Left boundary
1722+
# i = 0
1723+
# j=j_ind
1724+
# k=k_ind
1725+
# q = q[-1]+int_range(1,Ny*Nz)
1726+
# ii[q] = G[i,j,k].ravel()
1727+
# jj[q] = G[i,j,k].ravel()
1728+
# s[q] = 1.0
1729+
# q = q[-1]+int_range(1,Ny*Nz)
1730+
# ii[q] = G[i,j,k].ravel()
1731+
# jj[q] = G[i+1,j,k].ravel()
1732+
# s[q] = 1.0
1733+
# q = q[-1]+int_range(1,Ny*Nz)
1734+
# ii[q] = G[i,j,k].ravel()
1735+
# jj[q] = G[Nx,j,k].ravel()
1736+
# s[q] = -1.0
1737+
# q = q[-1]+int_range(1,Ny*Nz)
1738+
# ii[q] = G[i,j,k].ravel()
1739+
# jj[q] = G[Nx+1,j,k].ravel()
1740+
# s[q] = -1.0
1741+
# BCRHS[G[i,j,k].ravel()] = 0.0
1742+
#
1743+
#
1744+
#
17231745
if (not BC.front.periodic) and (not BC.back.periodic):
17241746
# Front boundary
17251747
k=Nz+1
@@ -1943,52 +1965,57 @@ def boundaryConditionsTermSpherical3D(BC: BoundaryConditions3D):
19431965
s[q] = -(BC.left.b/2 - BC.left.a/dx_1).ravel()
19441966
BCRHS[G[i,j,k].ravel()] = -(BC.left.c).ravel()
19451967
elif BC.right.periodic or BC.left.periodic: # periodic
1946-
# for a spherical coordinate system, the left and right boundaries (in the radial direction) cannot be periodic?
1947-
# or at least I cannot imagine them being periodic
1948-
# TODO: add a warning here; do the same for all radial boundaries
1949-
# Right boundary
1950-
i=Nx+1
1951-
j=j_ind
1952-
k=k_ind
1953-
q = q[-1]+int_range(1,Ny*Nz)
1954-
ii[q] = G[i,j,k].ravel()
1955-
jj[q] = G[i,j,k].ravel()
1956-
s[q] = 1.0
1957-
q = q[-1]+int_range(1,Ny*Nz)
1958-
ii[q] = G[i,j,k].ravel()
1959-
jj[q] = G[i-1,j,k].ravel()
1960-
s[q] = -1.0
1961-
q = q[-1]+int_range(1,Ny*Nz)
1962-
ii[q] = G[i,j,k].ravel()
1963-
jj[q] = G[0,j,k].ravel()
1964-
s[q] = dx_end/dx_1
1965-
q = q[-1]+int_range[1,Ny*Nz]
1966-
ii[q] = G[i,j,k].ravel()
1967-
jj[q] = G[1,j,k].ravel()
1968-
s[q] = -dx_end/dx_1
1969-
BCRHS[G[i,j,k].ravel()] = 0.0
1970-
1971-
# Left boundary
1972-
i = 0
1973-
j=j_ind
1974-
k=k_ind
1975-
q = q[-1]+int_range(1,Ny*Nz)
1976-
ii[q] = G[i,j,k].ravel()
1977-
jj[q] = G[i,j,k].ravel()
1978-
s[q] = 1.0
1979-
q = q[-1]+int_range(1,Ny*Nz)
1980-
ii[q] = G[i,j,k].ravel()
1981-
jj[q] = G[i+1,j,k].ravel()
1982-
s[q] = 1.0
1983-
q = q[-1]+int_range(1,Ny*Nz)
1984-
ii[q] = G[i,j,k].ravel()
1985-
jj[q] = G[Nx,j,k].ravel()
1986-
s[q] = -1.0
1987-
q = q[-1]+int_range(1,Ny*Nz)
1988-
ii[q] = G[i,j,k].ravel()
1989-
jj[q] = G[Nx+1,j,k].ravel()
1990-
s[q] = -1.0
1991-
BCRHS[G[i,j,k].ravel()] = 0.0
1968+
raise ValueError("Radial periodic boundary conditions are not physically meaningful.")
1969+
#
1970+
# Keep the following code for future reference, once a physically relevant
1971+
# case has been identified for radial periodic BCs...
1972+
#
1973+
# # Right boundary
1974+
# i=Nx+1
1975+
# j=j_ind
1976+
# k=k_ind
1977+
# q = q[-1]+int_range(1,Ny*Nz)
1978+
# ii[q] = G[i,j,k].ravel()
1979+
# jj[q] = G[i,j,k].ravel()
1980+
# s[q] = 1.0
1981+
# q = q[-1]+int_range(1,Ny*Nz)
1982+
# ii[q] = G[i,j,k].ravel()
1983+
# jj[q] = G[i-1,j,k].ravel()
1984+
# s[q] = -1.0
1985+
# q = q[-1]+int_range(1,Ny*Nz)
1986+
# ii[q] = G[i,j,k].ravel()
1987+
# jj[q] = G[0,j,k].ravel()
1988+
# s[q] = dx_end/dx_1
1989+
# q = q[-1]+int_range(1,Ny*Nz)
1990+
# ii[q] = G[i,j,k].ravel()
1991+
# jj[q] = G[1,j,k].ravel()
1992+
# s[q] = -dx_end/dx_1
1993+
# BCRHS[G[i,j,k].ravel()] = 0.0
1994+
1995+
# # Left boundary
1996+
# i = 0
1997+
# j=j_ind
1998+
# k=k_ind
1999+
# q = q[-1]+int_range(1,Ny*Nz)
2000+
# ii[q] = G[i,j,k].ravel()
2001+
# jj[q] = G[i,j,k].ravel()
2002+
# s[q] = 1.0
2003+
# q = q[-1]+int_range(1,Ny*Nz)
2004+
# ii[q] = G[i,j,k].ravel()
2005+
# jj[q] = G[i+1,j,k].ravel()
2006+
# s[q] = 1.0
2007+
# q = q[-1]+int_range(1,Ny*Nz)
2008+
# ii[q] = G[i,j,k].ravel()
2009+
# jj[q] = G[Nx,j,k].ravel()
2010+
# s[q] = -1.0
2011+
# q = q[-1]+int_range(1,Ny*Nz)
2012+
# ii[q] = G[i,j,k].ravel()
2013+
# jj[q] = G[Nx+1,j,k].ravel()
2014+
# s[q] = -1.0
2015+
# BCRHS[G[i,j,k].ravel()] = 0.0
2016+
#
2017+
#
2018+
#
19922019
if (not BC.front.periodic) and (not BC.back.periodic):
19932020
# Front boundary
19942021
k=Nz+1
@@ -2068,6 +2095,8 @@ def boundaryConditionsTermSpherical3D(BC: BoundaryConditions3D):
20682095
shape=((Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2)))
20692096
return BCMatrix, BCRHS
20702097

2098+
2099+
20712100
def boundaryConditionsTerm(BC):
20722101
"""
20732102
Generate the terms of the matrix equation representing the boundary conditions

0 commit comments

Comments
 (0)