1+ import numpy as np
2+ import matplotlib .pyplot as plt
3+ import treelog
4+ from nutils import mesh , function , solver , cli
5+ import precice
6+
7+ def main (nelems = 200 , dt = .005 , refdensity = 1e3 , refpressure = 101325.0 , psi = 1e-6 , viscosity = 1e-3 , theta = 0.5 ):
8+
9+ # --- preCICE initialization ---
10+
11+ participant_name = "Fluid1DLeft"
12+ config_file_name = "../precice-config.xml"
13+ solver_process_index = 0
14+ solver_process_size = 1
15+ participant = precice .Participant (participant_name , config_file_name , solver_process_index , solver_process_size )
16+ mesh_name = "Fluid1DLeft-Mesh"
17+ velocity_name = "Velocity"
18+ pressure_name = "Pressure"
19+ positions = [[0.0 , 500.0 , 0.0 ]]
20+ vertex_ids = participant .set_mesh_vertices (mesh_name , positions )
21+
22+ participant .initialize ()
23+ precice_dt = participant .get_max_time_step_size ()
24+
25+ # --- Nutils domain setup ---
26+
27+ domain , geom = mesh .rectilinear ([np .linspace (0 , 500 , nelems + 1 )]) # 1D domain from 0 to 500 with nelems elements
28+
29+ ns = function .Namespace ()
30+ ns .x = geom
31+ ns .dt = dt
32+ ns .θ = theta
33+ ns .ρref = refdensity
34+ ns .pref = refpressure
35+ ns .pin = 98100 # Inlet pressure (Pa)
36+ ns .μ = viscosity
37+ ns .ψ = psi
38+
39+ # Define basis functions: velocity (vector, degree 2), density (scalar, degree 1)
40+ ns .ubasis , ns .ρbasis = function .chain ([
41+ domain .basis ('std' , degree = 2 ).vector (domain .ndims ),
42+ domain .basis ('std' , degree = 1 )
43+ ])
44+
45+ # Define trial and previous-step functions
46+ ns .u_i = 'ubasis_ni ?lhs_n'
47+ ns .u0_i = 'ubasis_ni ?lhs0_n'
48+ ns .ρ = 'ρref + ρbasis_n ?lhs_n'
49+ ns .ρ0 = 'ρref + ρbasis_n ?lhs0_n'
50+ ns .p = 'pref + (ρ - ρref) / ψ'
51+ ns .utheta_i = 'θ u_i + (1 - θ) u0_i'
52+ ns .ρtheta = 'θ ρ + (1 - θ) ρ0'
53+
54+ # Cauchy stress tensor: viscous + pressure
55+ ns .σ_ij = 'μ (utheta_i,j + utheta_j,i) - p δ_ij'
56+
57+ # --- Weak form residuals ---
58+
59+ # Mass conservation
60+ res = domain .integral (
61+ 'ρbasis_n ((ρ - ρ0) / dt + (ρtheta utheta_k)_,k) d:x' @ ns ,
62+ degree = 4
63+ )
64+
65+ # Momentum conservation: unsteady + convective + diffusive terms
66+ res += domain .integral (
67+ '(ubasis_ni (((ρ u_i - ρ0 u0_i) / dt) + (ρtheta utheta_i utheta_j)_,j) + ubasis_ni,j σ_ij) d:x' @ ns ,
68+ degree = 4
69+ )
70+
71+ # Weakly enforced inlet pressure boundary condition
72+ res += domain .boundary ['left' ].integral (
73+ 'pin ubasis_ni n_i d:x' @ ns ,
74+ degree = 4
75+ )
76+
77+ # Initial condition
78+ lhs0 = np .zeros (res .shape )
79+
80+ # Time-stepping
81+ t = 0.0
82+ timestep = 0
83+ bezier = domain .sample ('bezier' , 2 )
84+
85+ f = open ("watchpoint.txt" , "w" )
86+
87+ # --- Time loop with preCICE coupling ---
88+ while participant .is_coupling_ongoing ():
89+
90+ # Read velocity data from other solver via preCICE
91+ u_read = participant .read_data (mesh_name , velocity_name , vertex_ids , precice_dt )
92+
93+ # Constrain outlet velocity to match coupled value
94+ stringintegral = f'(u_0 - ({ u_read [0 ][1 ]}
95+ sqr = domain .boundary ['right' ].integral (stringintegral @ ns , degree = 4 )
96+ cons = solver .optimize ('lhs' , sqr , droptol = 1e-14 )
97+
98+ # Write checkpoint if required by preCICE
99+ if participant .requires_writing_checkpoint ():
100+ lhscheckpoint = lhs0
101+
102+ # Solve nonlinear system (Newton's method)
103+ with treelog .context ('stokes' ):
104+ lhs = solver .newton (
105+ 'lhs' ,
106+ residual = res ,
107+ constrain = cons ,
108+ arguments = dict (lhs0 = lhs0 ),
109+ lhs0 = lhs0
110+ ).solve (1e-08 )
111+
112+ # Evaluate fields for visualization/debugging
113+ x , p , u , ρ = bezier .eval (['x_i' ,'p' , 'u_i' , 'ρ' ] @ ns , arguments = dict (lhs = lhs ))
114+
115+ # Send pressure at the right boundary to the other solver
116+ write_press = [[p [- 1 ]]]
117+ participant .write_data (mesh_name , pressure_name , vertex_ids , write_press )
118+
119+ # Advance in pseudo-time
120+ participant .advance (precice_dt )
121+ precice_dt = participant .get_max_time_step_size ()
122+
123+ # Restore checkpoint if needed
124+ if participant .requires_reading_checkpoint ():
125+ lhs0 = lhscheckpoint
126+ else :
127+ # Update old solution
128+ lhs0 = lhs
129+ timestep += timestep
130+
131+ # Save probe values (time, inlet pressure, inlet velocity, outlet pressure, outlet velocity, pressure at the middle, velocity at the middle)
132+ x , p , ρ , u = bezier .eval (['x_i' , 'p' , 'ρ' , 'u_i' ] @ ns , lhs = lhs )
133+ f .write ("%e; %e; %e; %e; %e; %e; %e\n " % (t , p [0 ], u [0 ], p [- 1 ], u [- 1 ], p [199 ], u [199 ]))
134+ f .flush ()
135+
136+ t += precice_dt # advance pseudo-time
137+
138+ # Finalize preCICE
139+ participant .finalize ()
140+ f .close ()
141+
142+ if __name__ == '__main__' :
143+ cli .run (main )
0 commit comments