Skip to content

Commit 1a777e2

Browse files
sbryngelsonsbryngelson
committed
reorganize and refactor
1 parent 5c686b6 commit 1a777e2

12 files changed

+737
-328
lines changed

‎simulate.m‎

Lines changed: 42 additions & 319 deletions
Large diffs are not rendered by default.

‎src/NS_2d_cylinder_PHS.m‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1-
function [W3,p] = NS_2d_cylinder_PHS(dt,nu,W1,W2,Dy,Dx,L_inv,L_u_inv,L_v_inv,L0,L_B,L_B_obs,L_W,L_B_y,L_B_S,D0_12_x,D0_12_y,D0_21_x,D0_21_y,Dy_b,Dy_b_1,D0_12_x_obs,D0_12_y_obs,p0,W0)
2-
%NS_2D_CYLINDER_PHS Fractional step method for 2D Navier-Stokes equations around cylinder
1+
function [W3,p] = NS_2d_fractional_step_PHS(dt,nu,W1,W2,Dy,Dx,L_inv,L_u_inv,L_v_inv,L0,L_B,L_B_obs,L_W,L_B_y,L_B_S,D0_12_x,D0_12_y,D0_21_x,D0_21_y,Dy_b,Dy_b_1,D0_12_x_obs,D0_12_y_obs,p0,W0)
2+
%NS_2D_FRACTIONAL_STEP_PHS Fractional step method for 2D incompressible Navier-Stokes equations
33
%
44
% This function implements the fractional step method for incompressible Navier-Stokes:
55
% 1. Advection-diffusion step using Adams-Bashforth + Crank-Nicolson
@@ -61,15 +61,15 @@
6161
V1 = V + dt * (ADAMS_BASHFORTH_COEFF_CURRENT*H_V - ADAMS_BASHFORTH_COEFF_PREVIOUS*H_V_1);
6262

6363
% Apply boundary conditions after advection step
64-
% Cylinder and inlet boundaries: maintain previous values (no-slip, prescribed inlet)
65-
U1(end-L_B+1:end-0) = U(end-L_B+1:end-0); % Cylinder + inlet BCs for u
64+
% Obstacle and inlet boundaries: maintain previous values (no-slip, prescribed inlet)
65+
U1(end-L_B+1:end-0) = U(end-L_B+1:end-0); % Obstacle + inlet BCs for u
6666

6767
% Wall boundaries: enforce du/dy = 0 (slip condition for u-velocity)
6868
U1(L_W+1:L_W+L_B_y) = -(Dy_b*U1)./Dy_b_1; % Wall BC: du/dy = 0
6969

7070
% v-velocity boundary conditions
7171
V1(L_W+1:L_W+L_B_y) = V(L_W+1:L_W+L_B_y); % Wall BC: v = 0 (no penetration)
72-
V1(end-L_B+1:end-0) = V(end-L_B+1:end-0); % Cylinder + inlet BCs for v
72+
V1(end-L_B+1:end-0) = V(end-L_B+1:end-0); % Obstacle + inlet BCs for v
7373
%
7474
%
7575
%% STEP 2: Viscous diffusion using Crank-Nicolson method
@@ -80,8 +80,8 @@
8080
V2 = V1 + dt*nu*(L0*V)*CRANK_NICOLSON_COEFF; % Explicit viscous term for v
8181

8282
% Apply boundary conditions for right-hand side
83-
U2(end-L_B+1:end) = U(end-L_B+1:end); % Cylinder + inlet BCs for u
84-
V2(end-L_B+1:end) = V(end-L_B+1:end); % Cylinder + inlet BCs for v
83+
U2(end-L_B+1:end) = U(end-L_B+1:end); % Obstacle + inlet BCs for u
84+
V2(end-L_B+1:end) = V(end-L_B+1:end); % Obstacle + inlet BCs for v
8585

8686
% Set boundary node values in RHS vector
8787
U2(L_W+1:end-L_B) = zeros(L_W2/2-L_W-L_B,1); % Outlet + wall BCs for u
@@ -127,12 +127,12 @@
127127
% Apply boundary conditions to final velocity field
128128

129129
% u-velocity boundary conditions
130-
U3(end-L_B+1:end) = U(end-L_B+1:end); % Cylinder + inlet BCs: u = prescribed
130+
U3(end-L_B+1:end) = U(end-L_B+1:end); % Obstacle + inlet BCs: u = prescribed
131131
U3(L_W+1:L_W+L_B_y) = -(Dy_b*U3)./Dy_b_1; % Wall BC: du/dy = 0
132132

133133
% v-velocity boundary conditions
134134
V3(L_W+1:L_W+L_B_y) = V(L_W+1:L_W+L_B_y); % Wall BC: v = 0 (no penetration)
135-
V3(end-L_B+1:end) = V(end-L_B+1:end); % Cylinder + inlet BCs: v = prescribed
135+
V3(end-L_B+1:end) = V(end-L_B+1:end); % Obstacle + inlet BCs: v = prescribed
136136

137137
%% Assemble final velocity vector for next time step
138138
W3 = [U3; V3]; % Combined velocity vector [u^(n+1); v^(n+1)]

‎src/NS_2d_fractional_step_PHS.m‎

Lines changed: 140 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
1+
function [W3,p] = NS_2d_fractional_step_PHS(dt,nu,W1,W2,Dy,Dx,L_inv,L_u_inv,L_v_inv,L0,L_B,L_B_obs,L_W,L_B_y,L_B_S,D0_12_x,D0_12_y,D0_21_x,D0_21_y,Dy_b,Dy_b_1,D0_12_x_obs,D0_12_y_obs,p0,W0)
2+
%NS_2D_FRACTIONAL_STEP_PHS Fractional step method for 2D incompressible Navier-Stokes equations
3+
%
4+
% This function implements the fractional step method for incompressible Navier-Stokes:
5+
% 1. Advection-diffusion step using Adams-Bashforth + Crank-Nicolson
6+
% 2. Pressure correction step to enforce incompressibility
7+
% 3. Velocity correction using pressure gradient
8+
%
9+
% INPUTS:
10+
% dt - Time step size
11+
% nu - Kinematic viscosity (1/Reynolds number)
12+
% W1, W2 - Velocity fields at previous time steps [U1; V1], [U2; V2]
13+
% Dy, Dx - Spatial derivative operators (y and x directions)
14+
% L_inv - Pressure Poisson solver (LU factorized)
15+
% L_u_inv, L_v_inv - Velocity diffusion solvers (LU factorized)
16+
% L0 - Laplacian operator for velocity diffusion
17+
% L_B, L_B_obs, L_W, L_B_y, L_B_S - Boundary indexing parameters
18+
% D0_12_x, D0_12_y - Divergence operators (V-grid to P-grid)
19+
% D0_21_x, D0_21_y - Gradient operators (P-grid to V-grid)
20+
% Dy_b, Dy_b_1 - Boundary derivative operators
21+
% D0_12_x_obs, D0_12_y_obs - Obstacle boundary operators
22+
% p0 - Previous pressure field
23+
% W0 - Initial velocity field (for reference)
24+
%
25+
% OUTPUTS:
26+
% W3 - Updated velocity field [U3; V3] at next time step
27+
% p - Updated pressure field
28+
29+
%% Load numerical scheme coefficients from configuration
30+
cfg = config();
31+
ADAMS_BASHFORTH_COEFF_CURRENT = cfg.schemes.adams_bashforth_current; % Coefficient for current time step
32+
ADAMS_BASHFORTH_COEFF_PREVIOUS = cfg.schemes.adams_bashforth_previous; % Coefficient for previous time step
33+
CRANK_NICOLSON_COEFF = cfg.schemes.crank_nicolson; % Coefficient for implicit diffusion
34+
35+
%% Extract velocity components from input vectors
36+
L_W2 = length(W2); % Total length of velocity vector (2 * number of nodes)
37+
38+
% Current time step velocities (time level n)
39+
U = W2(1:L_W2/2); % u-velocity component at current time
40+
V = W2(L_W2/2+1:end); % v-velocity component at current time
41+
42+
% Previous time step velocities (time level n-1)
43+
U_1 = W1(1:L_W2/2); % u-velocity component at previous time
44+
V_1 = W1(L_W2/2+1:end); % v-velocity component at previous time
45+
46+
47+
%% STEP 1: Advection using Adams-Bashforth method
48+
% Compute nonlinear advection terms: -u*du/dx - v*du/dy, -u*dv/dx - v*dv/dy
49+
50+
% Advection terms at current time step (time level n)
51+
H_U = (-U.*(Dx*(U)) - V.*(Dy*(U))); % Advection of u-momentum: -(u*du/dx + v*du/dy)
52+
H_V = (-U.*(Dx*(V)) - V.*(Dy*(V))); % Advection of v-momentum: -(u*dv/dx + v*dv/dy)
53+
54+
% Advection terms at previous time step (time level n-1)
55+
H_U_1 = (-U_1.*(Dx*(U_1)) - V_1.*(Dy*(U_1))); % Previous u-momentum advection
56+
H_V_1 = (-U_1.*(Dx*(V_1)) - V_1.*(Dy*(V_1))); % Previous v-momentum advection
57+
58+
% Adams-Bashforth extrapolation for advection terms
59+
% u^* = u^n + dt * (3/2 * H^n - 1/2 * H^(n-1))
60+
U1 = U + dt * (ADAMS_BASHFORTH_COEFF_CURRENT*H_U - ADAMS_BASHFORTH_COEFF_PREVIOUS*H_U_1);
61+
V1 = V + dt * (ADAMS_BASHFORTH_COEFF_CURRENT*H_V - ADAMS_BASHFORTH_COEFF_PREVIOUS*H_V_1);
62+
63+
% Apply boundary conditions after advection step
64+
% Obstacle and inlet boundaries: maintain previous values (no-slip, prescribed inlet)
65+
U1(end-L_B+1:end-0) = U(end-L_B+1:end-0); % Obstacle + inlet BCs for u
66+
67+
% Wall boundaries: enforce du/dy = 0 (slip condition for u-velocity)
68+
U1(L_W+1:L_W+L_B_y) = -(Dy_b*U1)./Dy_b_1; % Wall BC: du/dy = 0
69+
70+
% v-velocity boundary conditions
71+
V1(L_W+1:L_W+L_B_y) = V(L_W+1:L_W+L_B_y); % Wall BC: v = 0 (no penetration)
72+
V1(end-L_B+1:end-0) = V(end-L_B+1:end-0); % Obstacle + inlet BCs for v
73+
%
74+
%
75+
%% STEP 2: Viscous diffusion using Crank-Nicolson method
76+
% Treat diffusion terms implicitly for stability: u** = u* + dt*nu*del^2*u
77+
78+
% Add explicit part of viscous terms (from previous time step)
79+
U2 = U1 + dt*nu*(L0*U)*CRANK_NICOLSON_COEFF; % Explicit viscous term for u
80+
V2 = V1 + dt*nu*(L0*V)*CRANK_NICOLSON_COEFF; % Explicit viscous term for v
81+
82+
% Apply boundary conditions for right-hand side
83+
U2(end-L_B+1:end) = U(end-L_B+1:end); % Obstacle + inlet BCs for u
84+
V2(end-L_B+1:end) = V(end-L_B+1:end); % Obstacle + inlet BCs for v
85+
86+
% Set boundary node values in RHS vector
87+
U2(L_W+1:end-L_B) = zeros(L_W2/2-L_W-L_B,1); % Outlet + wall BCs for u
88+
V2(L_W+L_B_y+1:end-L_B) = zeros(L_W2/2-L_W-L_B-L_B_y,1); % Outlet BCs for v
89+
90+
U2(L_W+1:L_W+L_B_y) = zeros(L_B_y,1); % Wall BC for u
91+
V2(L_W+1:L_W+L_B_y) = zeros(L_B_y,1); % Wall BC for v
92+
93+
% Apply pressure-dependent boundary conditions on obstacle (from previous pressure)
94+
[L_B_obs_local,~] = size(D0_12_x_obs); % Number of obstacle boundary nodes
95+
U2(end-L_B_obs_local+1:end) = D0_12_x_obs*p0; % Obstacle BC with pressure: du/dn = dp/dx
96+
V2(end-L_B_obs_local+1:end) = D0_12_y_obs*p0; % Obstacle BC with pressure: dv/dn = dp/dy
97+
98+
% Solve implicit viscous step: (I - dt*nu/2*del^2) * u** = RHS
99+
U2 = L_u_inv(U2); % Solve for u-velocity after viscous diffusion
100+
V2 = L_v_inv(V2); % Solve for v-velocity after viscous diffusion
101+
%% STEP 3: Pressure correction to enforce incompressibility
102+
% Solve Poisson equation for pressure: del^2 p = div(u**)/dt
103+
104+
% Compute divergence of intermediate velocity field u**
105+
F0 = (D0_12_x*(U2-0) + D0_12_y*(V2-0)); % Divergence: du**/dx + dv**/dy
106+
107+
% Add boundary conditions and regularization to pressure system
108+
F = [F0; zeros(L_B_S+1,1)]; % Add zeros for boundary conditions and regularization
109+
110+
% Solve pressure Poisson equation: del^2 p = div(u**)/dt
111+
p = L_inv(F); % Solve using precomputed LU factorization
112+
113+
% Extract pressure field (remove regularization constraint)
114+
p = p(1:(length(F0)+L_B_S)); % Remove regularization component
115+
F = F(1:(length(F0)+L_B_S)); % Remove regularization component
116+
117+
118+
119+
%% STEP 4: Velocity correction using pressure gradient
120+
% Apply pressure correction: u^(n+1) = u** - dt * grad(p)
121+
122+
% Subtract pressure gradient from intermediate velocities
123+
U3 = (U2-0) - [(D0_21_x*p); zeros(L_B,1)]; % u^(n+1) = u** - dt*dp/dx
124+
V3 = (V2-0) - [(D0_21_y*p); zeros(L_B,1)]; % v^(n+1) = v** - dt*dp/dy
125+
126+
%% Final boundary condition enforcement
127+
% Apply boundary conditions to final velocity field
128+
129+
% u-velocity boundary conditions
130+
U3(end-L_B+1:end) = U(end-L_B+1:end); % Obstacle + inlet BCs: u = prescribed
131+
U3(L_W+1:L_W+L_B_y) = -(Dy_b*U3)./Dy_b_1; % Wall BC: du/dy = 0
132+
133+
% v-velocity boundary conditions
134+
V3(L_W+1:L_W+L_B_y) = V(L_W+1:L_W+L_B_y); % Wall BC: v = 0 (no penetration)
135+
V3(end-L_B+1:end) = V(end-L_B+1:end); % Obstacle + inlet BCs: v = prescribed
136+
137+
%% Assemble final velocity vector for next time step
138+
W3 = [U3; V3]; % Combined velocity vector [u^(n+1); v^(n+1)]
139+
140+
end

‎src/build_geometry.m‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
function G = build_geometry(cfg)
2+
%BUILD_GEOMETRY Generate geometry using appropriate geometry helper
3+
%
4+
% This function dispatches to the appropriate geometry generation function
5+
% based on the configuration and prints a summary of the selected geometry.
6+
%
7+
% INPUTS:
8+
% cfg - Configuration structure from config()
9+
%
10+
% OUTPUTS:
11+
% G - Geometry structure containing mesh, boundaries, and geometry-specific data
12+
13+
% Generate geometry using appropriate geometry helper
14+
% This supports different obstacle geometries via the config.geometry.type parameter
15+
switch lower(cfg.geometry.type)
16+
case 'cylinder'
17+
fprintf('Using cylinder geometry (radius=%.2f)...\n', cfg.geometry.obstacle_radius);
18+
G = make_cylinder_geometry(cfg);
19+
case 'ellipse'
20+
fprintf('Using ellipse geometry (a=%.2f, b=%.2f)...\n', cfg.geometry.ellipse_a, cfg.geometry.ellipse_b);
21+
G = make_ellipse_geometry(cfg);
22+
otherwise
23+
error('Unknown geometry type: %s. Supported types: cylinder, ellipse', cfg.geometry.type);
24+
end
25+
26+
end

‎src/build_intergrid_ops.m‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
function [D0_21_x, D0_21_y, D0_12_x, D0_12_y] = build_intergrid_ops(G, xy, xy1, xy_s, xy1_s, S, cfg)
2+
%BUILD_INTERGRID_OPS Generate interpolation operators between grids
3+
%
4+
% This function builds RBF-FD operators for interpolation and differentiation
5+
% between the velocity grid (V-grid) and pressure grid (P-grid).
6+
%
7+
% INPUTS:
8+
% G - Geometry structure
9+
% xy - Interior velocity nodes
10+
% xy1 - Complete velocity grid (interior + boundary)
11+
% xy_s - Interior pressure nodes
12+
% xy1_s - Complete pressure grid (interior + boundary)
13+
% S - Stencil structure from build_stencils
14+
% cfg - Configuration structure
15+
%
16+
% OUTPUTS:
17+
% D0_21_x - x-derivative operator (P-grid to V-grid)
18+
% D0_21_y - y-derivative operator (P-grid to V-grid)
19+
% D0_12_x - x-derivative operator (V-grid to P-grid)
20+
% D0_12_y - y-derivative operator (V-grid to P-grid)
21+
22+
% Differentiation operators: P-grid to V-grid
23+
% Low-order PHS-RBFs and polynomials near boundaries
24+
[D0_21_all] = RBF_PHS_FD_all(xy1(1:length(xy)+length(G.boundary_y)+length(G.boundary_out),:), xy1_s, S.Nearest_Idx_interp_21, cfg.rbf.order_interpolation_low, cfg.rbf.poly_degree_main, cfg.rbf.poly_degree_interpolation_low);
25+
26+
D0_21_x = D0_21_all{1};
27+
D0_21_y = D0_21_all{2};
28+
29+
% Use precomputed indices for nodes far from all boundaries (high-order region)
30+
Nearest_Idx_nc = find(G.idx_far_boundaries_V);
31+
xy_nc = xy1(Nearest_Idx_nc,:);
32+
33+
[D0_21_all_nc] = RBF_PHS_FD_all(xy_nc, xy1_s, S.Nearest_Idx_interp_21(Nearest_Idx_nc,:), cfg.rbf.order_interpolation_high, cfg.rbf.order_interpolation_high_poly, cfg.rbf.poly_degree_interpolation_high);
34+
35+
D0_21_x(Nearest_Idx_nc,:) = D0_21_all_nc{1};
36+
D0_21_y(Nearest_Idx_nc,:) = D0_21_all_nc{2};
37+
38+
% Differentiation operators: V-grid to P-grid
39+
% Low-order PHS-RBFs and polynomials near boundaries
40+
D0_12_all = RBF_PHS_FD_all(xy_s, xy1, S.Nearest_Idx_interp, cfg.rbf.order_boundary, cfg.rbf.order_interpolation_high_poly, cfg.rbf.derivative_order);
41+
42+
D0_12_x = D0_12_all{1};
43+
D0_12_y = D0_12_all{2};
44+
45+
% Use precomputed indices for pressure nodes far from all boundaries (high-order region)
46+
Nearest_Idx_nc = find(G.idx_far_boundaries_P);
47+
xy_nc = xy1_s(Nearest_Idx_nc,:);
48+
49+
D0_12_all_nc = RBF_PHS_FD_all(xy_nc, xy1, S.Nearest_Idx_interp(Nearest_Idx_nc,:), cfg.rbf.order_interpolation_high, cfg.rbf.order_interpolation_high_poly, cfg.rbf.poly_degree_interpolation_high);
50+
51+
D0_12_x(Nearest_Idx_nc,:) = D0_12_all_nc{1};
52+
D0_12_y(Nearest_Idx_nc,:) = D0_12_all_nc{2};
53+
54+
end

‎src/build_pressure_system.m‎

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
function P = build_pressure_system(G, xy_s, xy1_s, cfg)
2+
%BUILD_PRESSURE_SYSTEM Generate pressure system and boundary conditions
3+
%
4+
% This function builds the complete pressure Poisson system including:
5+
% - Laplacian operator for interior nodes
6+
% - Boundary condition operators for all boundaries
7+
% - Assembled system matrix with regularization
8+
% - Precomputed LU decomposition for efficient solving
9+
%
10+
% INPUTS:
11+
% G - Geometry structure
12+
% xy_s - Interior pressure nodes
13+
% xy1_s - Complete pressure grid (interior + boundary)
14+
% cfg - Configuration structure
15+
%
16+
% OUTPUTS:
17+
% P - Structure containing pressure system components:
18+
% .L_inv_s - Pressure solver function handle
19+
% .D0_21_x_obs - Obstacle x-gradient operator
20+
% .D0_21_y_obs - Obstacle y-gradient operator
21+
% .L_s - Laplacian operator
22+
23+
% Extract boundary data from geometry
24+
boundary_y_s = G.boundary_y_s;
25+
boundary_out_s = G.boundary_out_s;
26+
boundary_in_s = G.boundary_in_s;
27+
boundary_obs_s = G.boundary_obs_s;
28+
boundary_s = [boundary_y_s; boundary_out_s; boundary_in_s; boundary_obs_s];
29+
30+
% Build k-nearest neighbor stencils for pressure grid
31+
k = cfg.rbf.stencil_size_main;
32+
Nearest_Idx_s = nearest_interp(xy1_s, xy1_s, k);
33+
34+
% Generate Laplacian operator for pressure Poisson equation (P-grid to P-grid)
35+
[D_s_all] = RBF_PHS_FD_all(xy_s, xy1_s, Nearest_Idx_s(1:length(xy_s),:), cfg.rbf.order_main, cfg.rbf.poly_degree_main, cfg.rbf.laplacian_order);
36+
P.L_s = D_s_all{3}; % Extract Laplacian operator (del^2 p = d^2p/dx^2 + d^2p/dy^2)
37+
38+
% Generate boundary condition operators for pressure system
39+
% Obstacle boundary: Neumann BC (dp/dn = 0, normal derivative = 0)
40+
[Dn1_b_s, Nearest_Idx_b_obs] = build_obstacle_pressure_bc(G, xy_s, xy1_s, cfg);
41+
42+
% Wall boundaries (top/bottom): Neumann BC (dp/dy = 0)
43+
[Nearest_Idx_b_y] = nearest_interp(boundary_y_s, [xy_s; xy1_s(length(xy_s)+length(boundary_y_s)+1,:)], cfg.rbf.stencil_size_boundary_wall);
44+
Nearest_Idx_b_y = [(length(xy_s)+1:length(xy_s)+length(boundary_y_s))' Nearest_Idx_b_y];
45+
46+
D = RBF_PHS_FD_all(boundary_y_s, xy1_s, Nearest_Idx_b_y, cfg.rbf.order_boundary, cfg.rbf.poly_degree_boundary, cfg.rbf.derivative_order);
47+
Dy_b_s = D{2}; % y-derivative operator for wall boundary condition
48+
49+
% Inlet boundary: Neumann BC (dp/dx = 0)
50+
[Nearest_Idx_b_x] = nearest_interp(boundary_in_s, xy1_s(1:length(xy_s)+length(boundary_y_s),:), cfg.rbf.stencil_size_boundary_wall);
51+
Nearest_Idx_b_x = [(length(xy1_s)+1-length(boundary_obs_s)-length(boundary_in_s):length(xy1_s)-length(boundary_obs_s))' Nearest_Idx_b_x];
52+
53+
D = RBF_PHS_FD_all(boundary_in_s, xy1_s, Nearest_Idx_b_x, cfg.rbf.order_boundary, cfg.rbf.poly_degree_boundary, cfg.rbf.derivative_order);
54+
Dx_in_s = D{1}; % x-derivative operator for inlet boundary condition
55+
56+
% Outlet boundary: Neumann BC (dp/dx = 0)
57+
[Nearest_Idx_b_x] = nearest_interp(boundary_out_s, xy1_s(1:length(xy_s)+length(boundary_y_s),:), cfg.rbf.stencil_size_boundary_outlet);
58+
Nearest_Idx_b_x = [(length(xy_s)+1+length(boundary_y_s):length(xy_s)+length(boundary_y_s)+length(boundary_out_s))' Nearest_Idx_b_x];
59+
60+
D = RBF_PHS_FD_all(boundary_out_s, xy1_s, Nearest_Idx_b_x, cfg.rbf.order_boundary, cfg.rbf.poly_degree_boundary, cfg.rbf.derivative_order);
61+
Dx_out_s = D{1}; % x-derivative operator for outlet boundary condition
62+
63+
% Assemble pressure boundary condition matrices
64+
% Initialize boundary condition matrices (each row corresponds to one boundary node)
65+
L_bc_c = zeros(length(boundary_s), length(xy1_s)); % Obstacle BC matrix
66+
L_bc_out = zeros(length(boundary_s), length(xy1_s)); % Outlet BC matrix
67+
L_bc_y = zeros(length(boundary_s), length(xy1_s)); % Wall BC matrix
68+
L_bc_in = zeros(length(boundary_s), length(xy1_s)); % Inlet BC matrix
69+
70+
% Fill obstacle boundary condition matrix (dp/dn = 0)
71+
Idx_boundary_obs = length(xy1_s)-length(boundary_obs_s)+[1:length(boundary_obs_s)];
72+
L_bc_c(Idx_boundary_obs-length(xy_s),:) = Dn1_b_s;
73+
74+
% Fill outlet boundary condition matrix (dp/dx = 0)
75+
Idx_boundary_out = length(xy_s) + length(boundary_y_s) + [1:length(boundary_out_s)];
76+
L_bc_out(Idx_boundary_out-length(xy_s),:) = Dx_out_s;
77+
78+
% Fill wall boundary condition matrix (dp/dy = 0)
79+
Idx_boundary_y = length(xy_s) + [1:length(boundary_y_s)];
80+
L_bc_y(Idx_boundary_y-length(xy_s),:) = Dy_b_s;
81+
82+
% Fill inlet boundary condition matrix (dp/dx = 0)
83+
Idx_boundary_in = length(xy1_s) - length(boundary_obs_s) - length(boundary_in_s) + [1:length(boundary_in_s)];
84+
L_bc_in(Idx_boundary_in-length(xy_s),:) = Dx_in_s;
85+
86+
% Assemble complete pressure system: [Laplacian at interior; BCs at boundary]
87+
L1 = [P.L_s(1:length(xy_s),:); L_bc_c+L_bc_out+L_bc_y+L_bc_in];
88+
89+
% Add regularization to fix pressure datum (pressure is defined up to constant)
90+
% This adds the constraint sum(p) = 0 to make the system uniquely solvable
91+
L1 = [L1 [ones(length(xy_s),1); zeros(length(boundary_s),1)]; [ones(1,length(xy_s)) zeros(1,length(boundary_s))] 0];
92+
93+
% Precompute LU decomposition for efficient pressure solves
94+
[LL,UU,pp,qq,rr] = lu(L1);
95+
P.L_inv_s = @(v) (qq*(UU\(LL\(pp*(rr\(v)))))); % Pressure solver function
96+
97+
% Generate obstacle boundary gradient operators for pressure-dependent boundary conditions
98+
[P.D0_21_x_obs, P.D0_21_y_obs] = build_obstacle_grad_operators(G, xy1_s, cfg);
99+
100+
end

0 commit comments

Comments
 (0)