Physics-integrated neural differentiable modeling for immersed boundary systems

Ph ysics-integrated neural differentiable modeling f or immersed boundar y sys tems Chenglin Li a , Hang Xu a , ∗ , Jianting Chen b and Y anf ei Zhang b a State Key Lab of Ocean Engineering, School of Ocean and Civil Engineer ing, Shanghai Jiao T ong Univer sity, Shanghai 200240, China b State Key Laborat or y of Maritime T ec hnology and Safety , Shanghai Ship and Shipping Resear ch Institute Co., Ltd, Shanghai 200135, China. A R T I C L E I N F O Keyw ords : Phy sical hybrid neural modeling Differentiable programming Immersed boundar y Long-term rolling for ecast A B S T R A C T Accuratel y , efficiently , and stably computing complex fluid flow s and their ev olution near solid boundaries ov er long horizons remains challenging. Conv entional numerical sol vers require fine g rids and small time steps to resolv e near-w all dynamics, resulting in high computational costs, while purely data-driven surrogate models accumulate rollout errors and lack robustness under extrapolativ e conditions. To address these issues, this study extends existing neural PDE solv ers by developing a ph ysics-integrated differentiable framew ork f or long-horizon prediction of immersed-boundary flow s. A key design aspect of the framework includes an important improv ement, namely the s tr uctural integ ration of physical pr inciples into an end-to-end dif- f erentiable architecture incorporating a PDE-based inter mediate velocity module and a multi- direct forcing immersed boundar y module, both adhering to the pressure-projection procedure f or incompressible flo w computation. The computationally e xpensive pressure projection step is substituted with a lear ned implicit cor rection using ConvR esNet blocks to reduce cost, and a sub- iteration strategy is introduced to separate the embedded physics module’ s stability requirement from t he sur rogate model’s time step, enabling stable coarse-gr id autoregressive rollouts with large effective time increments. The framewor k uses only single-step supervision f or training, eliminating long-hor izon backpropagation and reducing training time to under one hour on a single GPU. Evaluations on benchmark cases of flow past a stationary cylinder and a rotationally oscillating cylinder at 𝑅𝑒 = 100 show the proposed model consistentl y outperforms purel y data- driven, physics-loss-constrained, and coarse-gr id numer ical baselines in flow -field fidelity and long-horizon stability , while achie ving an appro ximately 200-fold infer ence speedup ov er the high-resolution solv er . 1. Introduction Complex fluid flow s and their ev olution in the presence of solid boundar ies represent a fundament al class of problems in engineering and natural sy stems. Despite significant advances in bo th accuracy and efficiency , con ventional computational fluid dynamics (CFD) still suffers fr om high computational costs, especiall y when fine resolutions of geometric boundaries and structural dynamic responses are required. This cost becomes par ticularl y prohibitive in flo w control and optimization [ 1 ], where large numbers of long-horizon simulations are typically requir ed to identify optimal parameters [ 2 ] and control strategies [ 3 , 4 ]. Hence, dev eloping efficient computational framew orks is of significant import ance and has broad practical relev ance. For flo w computations inv olving comple x and/or moving boundaries, the immersed boundar y (IB) method [ 5 ] solv es Euler ian field variables on a fixed Cartesian g rid and enforces boundary conditions by cor recting near -wall flow s using boundary/object inf or mation described in a Lag rangian manner . In this w ay , it av oids the repeated mesh regeneration associated with body-fitted [ 6 ] gr ids under boundary motion[ 7 ]. The ke y advantages of IB methods include highly regular computational procedures and data str uctures [ 8 , 9 ], fa vorable parallel scalability [ 10 , 11 ]. These properties enable high-throughput time adv ancement and f acilitate embedding boundary- handling modules as operators wit hin an end-to-end comput ational framew ork. With the dev elopment of scientific mac hine learning, neural netw orks ha ve been widel y used to construct f ast surrogate models for fluid [ 12 ] and fluid–structure sys tems [ 13 ]. Representativ e approaches include lear ning spatiotemporal ev olution in a latent space via proper or thogonal decomposition (POD) [ 14 ] or encoder–decoder architectures[ 15 ], and lear ning field mappings on structured/unstr uctured meshes using conv olutional neural networks (CNNs) [ 16 , 17 ] or graph neural netw orks (GNNs) [ 18 , 19 ]. These methods demonstrate t he potential of data-dr iv en prediction f or high-dimensional spatiotemporal dynamics. Ho wev er, man y purel y data-dr iv en black -bo x models rel y ∗ Corresponding author . E-mail address: hangxu@sjtu.edu.cn (H. Xu). OR CID (s): 0009-0004-3842-4391 (C. Li); 0000-0003-4176-0738 (H. Xu) CL Li, H Xu: Preprint submitted to Elsevier Page 1 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws hea vily on large datasets, and t hey of ten suffer from limited generalization under e xtrapolation in parameter space as well as er ror accumulation in long-hor izon rollouts [ 20 ]. Moreover , the cost of data acq uisition itself ma y become a major bottleneck for practical engineer ing deployment [ 21 ]. This has motiv ated a renew ed f ocus on lever aging accumulated phy sical priors (governing equations and con- straints) to reduce t he training difficulty of purely dat a-driven models [ 22 ], while exploiting the nonlinear representation capacity of neural netw orks to alleviate the computational burden of ph ysics-based solvers [ 23 ]. Ac hieving t his goal requires a careful balance between traditional physics-based framew orks and dat a-driven sur rogates. Existing physics- integrated architectures can be broadly categorized into two lines: (i) incor porating gov erning equations into the loss functions of deep neural networks (i.e., physics-inf or med neural ne tworks, PINNs) [ 22 ], and (ii) embedding ph ysical structure into neural architectures [ 24 ] (i.e., PDE-preser ved neural netw orks, ref er red to as PPNNs [ 25 ]). PINNs incorporate phy sical priors directl y into the training objective, typicall y by regular izing the netw ork with residuals of t he go verning eq uations to impro ve ph ysical consistency and mitigate t he lac k of constraints in purely data-dr iv en models. Since their introduction, PINNs ha ve been successfully applied across fluid dynamics [ 26 , 27 ], solid mechanics [ 28 ], heat transf er [ 29 ], and related areas. In fluid dynamics, PINNs ha ve been used f or sol ving the Na vier–Stokes equations, pro viding an alternative route f or rapid simulation of complex flow s. Ne vertheless, PINNs usually constr uct par tial differential equation (PDE) residuals via automatic differentiation (AD) [ 30 ] and/or numer ical- discretization-based reconstr uction, and include them as soft constraints in the loss. This often increases optimization difficulty and mak es training highl y sensitiv e to h yper parameters such as the relative weights of loss ter ms. In high- dimensional parameter spaces or under complex boundary conditions, balancing equation losses and dat a losses can furt her lead to training instability or limited gains [ 31 ]. The second research line, of ten ref er red to as PPNNs, embeds physical priors directly into the networ k architecture through discrete operators. By structurally constraining the hypothesis space, this design reduces the burden on purely data-dr iv en lear ning and mitigates t he accumulation and propagation of nonphysical er rors dur ing long-hor izon rollouts. The central idea is to establish a mathematical connection be tween neural-ne twor k structure and numerical PDE solvers [ 25 , 32 , 33 ]. This connection, together wit h differentiable programming [ 34 , 35 ], enables physical pr iors to be incorporated into the architecture itself. Such designs provide deeper insights into integrating phy sics pr iors wit h data-dr iv en surrogates. Related w orks hav e demonstrated long-time stable computations f or complex fluid problems [ 36 , 37 , 38 ], and ha ve embedded str uctural-dynamics solv ers to enable accurate prediction of structural responses and object boundar ies [ 36 ]. Although hybrid differentiable neural modeling has shown strong potential, the field remains at an early stage and requires fur ther dev elopment, particularly in reducing computational and inf erence costs for long-horizon prediction [ 25 , 36 , 39 ]. Based on the af orementioned literature sur ve y , long-hor izon prediction f or complex immersed boundary flo w problems still confronts two core challeng es that hinder its practical engineer ing application, and t hese challeng es motiv ate t he innov ations proposed in this w ork. Firs t, stable long-time f orecasts in existing machine learning-based approaches typically rely on sequence-modeling framewor ks and multi-step training strategies, which inevitably lead to a significant sur ge in training cost and time consumption [ 39 , 40 ]. Moreov er, phy sics-hybrid models, despite their advantages in er ror mitigation, often suffer from instability at the earl y training stage and ev en immediate diverg ence, posing a critical barr ier to their deployment [ 36 ]. Second, embedding physical structure into netw ork architectures, while constrains the hypo thesis space and reduces effectiv e model capacity , introduces inherent numerical restrictions inherited from traditional physics solv ers. If these restrictions are not properly addressed, phy sical constraints may tur n counterproductive b y undermining long-horizon stability , slowing do wn inference speed, and reducing deplo yment flexibility . This directl y conflicts with the core goal of dev eloping efficient computational frame wor ks f or comple x flow simulations. T o address t he af orementioned challenges, we extend the neural differentiable model introduced b y Fan and W ang [ 36 ] for long-horizon prediction of flow s subjected to immersed boundar y sy stems. Our approach is built upon the PPNNs paradigm [ 25 ], which establishes an e xplicit cor respondence betw een the netw ork architecture and numerical PDE discretizations. The work presents f our ke y innovations, each specifically t arg eting the limitations identified in existing research. For a star t, t he model architecture and dat a flow str ictly follo w the pressure-projection procedure f or incompressible flow s, ensur ing each intermediate variable has a clear physical meaning and thereby significantly mitigating the issue of nonph ysical error accumulation [ 20 ]. In addition, a sub-iteration strategy decouples the embedded phy sics solver ’ s s tability constraint from the model time step, which in turn enables stable coarse-grid rollouts wit h larg er time steps to reduce comput ational over head [ 36 ]. Another key innov ation lies in the replacement of t he traditional pressure-Poisson projection with a learned correction module, which retains incompressibility CL Li, H Xu: Preprint submitted to Elsevier P age 2 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws constraints while furt her reducing computational burden [ 22 , 21 ]. Lastl y , the model achie ves long-horizon st able prediction via single-step supervised training, whic h a voids long-hor izon bac kpropagation to effectivel y reduce both computational ov erhead and training complexity [ 39 , 40 ]. 2. Case setting T o evaluate the perf or mance of the proposed ph ysics-integrated neural model, w e consider tw o represent ativ e IB method benchmark cases: (i) flow past a stationar y circular cylinder , and (ii) flo w past a circular cy linder subjected to a prescr ibed rotational oscillation. The computational domain and boundary conditions for both cases are shown in Fig. 1 . Unless otherwise stated, the Reynolds number is fixed at 𝑅𝑒 = 100 , wit h an inflow v elocity of 𝑈 ∞ = 1 m/s and a cylinder diameter of 𝐷 = 1 m. First, high-fidelity simulations are conducted on a high-resolution gr id using a multi-step direct-forcing IBM scheme; t he cor responding numer ical v er ification is pro vided in the Appendix A . These simulations ser ve to constr uct the high-resolution dat aset used as ground tr uth f or both training and testing. In addition, a low -resolution training dataset is obtained b y do wnsampling the high-resolution solutions, with the do wnsampling procedure detailed in the Appendix B.1 . It is worth noting that, at the classical benchmark condition of 𝑅𝑒 = 100 , the flow past a stationary cylinder exhibits a pronounced per iodic vortex-shedding behavior , which to some extent reduces the difficulty of learning the unsteady flow ev olution and the associated hydrodynamic f orce responses. T o provide a more c hallenging e valuation scenario, we further consider the flo w past a cylinder subjected to rotational oscillation (see Fig. 1 ), where t he prescribed angular velocity is defined as f ollow s: 𝜔 𝑘 = 𝜔 𝑎 sin  2 𝜋 𝑓 𝑟 𝑘  , (1) where 𝜔 𝑎 denotes the angular -velocity amplitude, 𝑓 𝑟 is the prescribed (dimensionless) rotation frequency , and 𝑘 is t he 𝑘 -th learning time step. In the present w ork, we se t 𝑓 𝑟 = 1∕5 , i.e., the imposed rotational oscillation has a per iod of five lear ning time steps. Under rotational perturbations, the wake exhibits more intricate spatiotem poral e v olution. Mean while, the structural load response comprises not only components at the natural frequency but also modulation effects and additional spectral content associated with the imposed rotational-oscillation frequency [ 41 ], thereby markedly increasing the complexity of the data distribution and the difficulty of lear ning. Consequentl y , this benchmar k provides a more stringent test of the model’ s capability to represent unsteady flo ws and to generalize across complex dynamical regimes. 3. Methodology 3.1. Problem formulation The flow around the cylinder is gov er ned by the incompressible Na vier–Stokes (NS) equations: ∇ ⋅ 𝐮 = 0 , ( 𝐱 , 𝑡 ) ∈ Ω 𝑓 × [0 , 𝑇 ]; (2) 𝜕 𝐮 𝜕 𝑡 = −( 𝐮 ⋅ ∇) 𝐮 + 𝜈 ∇ 2 𝐮 − 1 𝜌 ∇ 𝑝, ( 𝐱 , 𝑡 ) ∈ Ω 𝑓 × [0 , 𝑇 ] . (3) Here 𝑡 and 𝐱 denote time and Euler ian spatial coordinates, respectiv ely . The velocity field 𝐮 ( 𝑡, 𝐱 ) and pressure 𝑝 ( 𝑡, 𝐱 ) are the primary spatiotemporal v ariables in Ω 𝑓 ⊂ ℝ 2 , while 𝜌 and 𝜈 are t he fluid density and kinematic viscosity . Given appropriate initial and boundar y conditions (IC/BCs), the velocity–pressure solution is well posed and therefore uniquel y determined. T o resol ve t he fluid–solid boundary interaction in this flow system, a multi-direct forcing immersed boundary method (IBM) is emplo yed, where the coupling between the Euler ian (fluid) and Lagrangian (solid) v ariables is computed via a Dirac delt a function. The det ailed coupling relations are given as follo ws: 𝜙 ( 𝐱 , 𝑡 ) = ∫ 𝑟 𝐱 Φ( 𝐗 , 𝑡 ) 𝛿 ( 𝐱 − 𝐗 ) 𝑑 𝛾 𝐗 , 𝛾 𝐗 ∈ Γ 𝐼 𝐵 ; (4) Φ( 𝐗 , 𝑡 ) = ∫ 𝑟 𝐱 𝜙 ( 𝐱 , 𝑡 ) 𝛿 ( 𝐱 − 𝐗 ) 𝑑 𝛾 𝐱 , 𝛾 𝐱 ∈ Γ 𝐼 𝐵 . (5) CL Li, H Xu: Preprint submitted to Elsevier P age 3 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Here Φ and 𝜙 are defined on the Lagrangian and Eulerian framew orks, respectively ; 𝐱 and 𝐗 cor respond to the se t of Eulerian grid nodes and Lagrangian mark ers in the immersed inter face boundary Γ , respectiv ely , and 𝛿 deno tes t he Dirac delta inter polation function. In our analysis, a three-point discrete formulation is adopted for this delt a function, as given by: 𝜙 ( 𝑟 ) =        1 3  1 +  1 − 3 𝑟 2  , 0 ≤ 𝑟 < 0 . 5 , 1 6  5 − 3 𝑟 −  −2 + 6 𝑟 − 3 𝑟 2  , 0 . 5 ≤ 𝑟 < 1 . 5 , 0 , 𝑟 ≥ 1 . 5 , (6) where 𝑟 =  𝑑  ∕ ℎ is the dimensionless normalized distance, 𝑑 denotes the ph ysical distance betw een a Lagrangian marker point and an Euler ian g rid node along one coordinate direction, and ℎ is the gr id spacing. 3.2. Ph ysics-integrated differentiable neural model The proposed IBM-based physics-integrated neural model, as shown in Fig. 2 , is constr ucted by f ollowing the dataflow of a classical pressure-projection sol ver , with the pressure-cor rection step perf ormed implicitly b y the netw ork module—a k ey distinction from traditional sol vers. The prediction of the intermediate velocity is built upon a principled connection betw een the neural architecture and the numer ical structure of the PDE: as discussed in [ 25 ], t he spatial discretization of differential operators can be inter preted as con volution wit h fix ed kernels. To ensure t he numerical stability of the embedded solv er, we further adopt a splitting strategy , in which t he intermediate-velocity update is car ried out via sub-iterations. Details of the fluid and str uctural solvers are pro vided as follo ws: The ov erall architecture of the model is primar ily organized to mirror the workflo w of a classical pressure- projection method. Specificall y , spatial discretization is encoded via finite differences on a staggered gr id, and temporal advancement is per formed using a first-order f or war d Euler scheme. To initiate the velocity prediction process, the intermediate velocity  𝐮 ∗ is first computed using the follo wing discrete formulation:  𝐮 ∗ = 𝐮 𝑡 + Δ 𝑡  −( 𝐮 𝑡 ⋅ ∇) 𝐮 𝑡 + 𝜈 ∇ 2 𝐮 𝑡  , (7) where Δ 𝑡 is the lear ning time step, which represents the time interval betw een tw o consecutive model predictions and is determined by t he temporal sampling inter v al of the training dataset. Notabl y , t he model’ s Δ 𝑡 is much larger than the underlying simulation time step 𝑑 𝑡 , which far ex ceeds the CFL stability limit, leading to unreliable solution of the phy sics-based constraint. This issue results in substantial numer ical er rors and spur ious noise, which sever ely deg rade the accuracy and stability of both training and inf erence. T o address this c hallenge and enable large-s tep inf erence while maintaining numerical robustness, the intermediate-velocity update is ref or mulated using a splitting s trategy . In practice, t his strategy is implemented through sub-iterations, with the iteration count 𝑁 determined according to t he CFL cons traint. For this work, 𝑁 = 20 , and according ly , the update in Eq. ( 7 ) is rewritten as: 𝐮 𝑛 +1 = 𝐮 𝑛 + 𝛿 𝑡  −( 𝐮 𝑛 ⋅ ∇) 𝐮 𝑛 + 𝜈 ∇ 2 𝐮 𝑛  , 𝑛 = 0 , 1 , … , 𝑁 − 1 , (8) where 𝛿 𝑡 = Δ 𝑡 ∕ 𝑁 is t he sub-iteration time s tep. In this work, Δ 𝑡 = 0 . 5 is determined b y the temporal sampling interval of t he training dataset, and 𝑁 = 20 sub-iterations are emplo yed such that 𝛿 𝑡 = Δ 𝑡 ∕ 𝑁 = 0 . 025 satisfies the CFL constraint on the coarse grid. When 𝑛 = 0 , the sub-iteration is initialized with the v elocity at the cur rent time le vel, i.e., 𝐮 0 = 𝐮 𝑡 ; after completing 𝑁 sub-steps, the inter mediate velocity is obtained as  𝐮 ∗ = 𝐮 𝑁 . F or spatial discretization, the con vectiv e ter ms are approximated using a first-order upwind scheme, while the diffusive ter ms are discretized wit h a second-order central difference scheme. The differential operators are applied using fixed-parameter con volution kernels, ensur ing st able and efficient updates for the velocity field during sub-iterations. This design not only preser v es the mat hematical st ability and physical fidelity of the imposed constraints but also prev ents the ov erall model from being bottlenecked by the numerical stability limit of the physics solver . Under this scheme,  𝐮 ∗ is computed b y the physics module on a sparse g rid via time-step sub-iterations, but it ma y still contain non-negligible er rors and numer ical noise. T o eliminate these residuals, a trainable ConvR esNet module is applied to cor rect  𝐮 ∗ , yielding the final inter mediate velocity  𝐮 :  𝐮 =  𝑐 𝑜𝑛𝑣 1 [  𝐮 ∗ , 𝐮 𝑡 ; 𝜃 1 ] , (9) CL Li, H Xu: Preprint submitted to Elsevier P age 4 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws where 𝜃 1 is the trainable parameters of t he ConvResN et module  𝑐 𝑜𝑛𝑣 1 . At this stage, the computed final intermediate v elocity  𝐮 has not ye t incor porated the influence of the immersed solid boundar ies wit hin the flo w domain. W e theref ore emplo y a multi-step direct-for cing immersed boundar y (IB) procedure to enforce the no-slip constraint on t he immersed inter face. Notably , all IB operations are per f ormed on the super-resol ved high-resolution fields (t he purple modules in Fig. 2 ). At the 𝑚 -t h IB sub-iteration, the Eulerian intermediate velocity  𝐮 𝑚 is first inter polated to t he Lag rangian marker locations 𝐗 𝐾 ∈ 𝛾 𝐗 through the regularized discrete delt a ker nel, yielding t he marker velocity :  𝐔 𝑚 ( 𝐗 𝐾 ) =  𝐱 𝑖 ∈ 𝛾 𝐱  𝐮 𝑚 ( 𝐱 𝑖 ) 𝛿 ℎ ( 𝐱 𝑖 − 𝐗 𝐾 ) ℎ 2 , 𝐗 𝐾 ∈ 𝛾 𝐗 , 𝑚 = 0 , 1 , … , 𝑀 − 1 , (10) where  𝐔 𝑚 ( 𝐗 𝐾 ) denotes t he inter polated velocity at the 𝐾 -th mark er, ℎ is the unif orm grid spacing, and 𝛿 ℎ ( ⋅ ) is the regular ized discrete delta k er nel. This kernel is constructed as the two-dimensional tensor product of t he one- dimensional kernel 𝜙 ( 𝑟 ) defined in Eq. ( 6 ), specifically: 𝛿 ℎ ( 𝐱 𝑖 − 𝐗 𝐾 ) = 1 ℎ 2 𝜙        𝑥 (1) 𝑖 − 𝑋 (1) 𝐾    ℎ     𝜙        𝑥 (2) 𝑖 − 𝑋 (2) 𝐾    ℎ     , (11) where 𝑥 (1) 𝑖 and 𝑥 (2) 𝑖 denote the two Cartesian components of the Eulerian node 𝐱 𝑖 , and 𝑋 (1) 𝐾 and 𝑋 (2) 𝐾 represent the two Cartesian components of the Lagrangian marker 𝐗 𝐾 . The factor 1∕ ℎ 2 ensures that the discre te kernel approximates the two-dimensional Dirac delt a distr ibution in the continuum limit. Based on the mismatch between the interpolated velocity  𝐔 𝑚 𝐾 and t he prescr ibed target boundary velocity 𝐔 𝑡𝑎𝑟𝑔 𝑒𝑡 𝐾 , t he direct-for cing term on each marker is defined as: 𝐅 𝑚 𝐾 = 𝐔 𝑡𝑎𝑟𝑔 𝑒𝑡 𝐾 −  𝐔 𝑚 𝐾 Δ 𝑡 , 𝐗 𝐾 ∈ 𝛾 𝐗 , 𝑚 = 0 , 1 , … , 𝑀 − 1 , (12) where Δ 𝑡 is the learning time step (i.e., the model’ s prediction inter val), which is consistent with the global step defined in Eq. ( 7 ). This consistency is enabled by our sub-iteration design, which decouples the ov erall model from the restr ictiv e stability constraint associated with a small simulation time step. After determining t he marker f orcing, it is subsequentl y spread back to the Euler ian g rid to form a body-f orce density: 𝐟 𝑚 ( 𝐱 𝑖 ) =  𝐗 𝐾 ∈ 𝛾 𝐗 𝐅 𝑚 𝐾 𝛿 ℎ ( 𝐱 𝑖 − 𝐗 𝐾 )Δ 𝑠 𝐾 , 𝐱 𝑖 ∈ 𝛾 𝐱 , (13) with Δ 𝑠 𝐾 = ℎ 𝛿𝑠 denoting t he discrete boundary measure associated wit h mark er 𝐾 , where 𝛿𝑠 is the mean spacing betw een adjacent Lag rangian markers. Using this body-f orce density , the Eulerian v elocity field is then corrected accordingly :  𝐮 𝑚 +1 ( 𝐱 𝑖 ) =  𝐮 𝑚 ( 𝐱 𝑖 ) + Δ 𝑡 𝐟 𝑚 ( 𝐱 𝑖 ) , 𝑚 = 0 , 1 , … , 𝑀 − 1 , (14) The abo ve gather-f orce-spread-update cycle (from Eq. 10 to Eq. 14 ) is repeated f or 𝑀 sub-iterations to ensure t he no-slip constraint is fully enforced. In t his work, 𝑀 = 5 sub-iterations are employ ed, which is f ound to be sufficient f or the no-slip constraint to conv erge on t he immersed boundary . Upon completion of the IB sub-iterations, the boundar y- cor rected final inter mediate velocity is obtained as: 𝐮 ∗ =  𝐮 𝑀 , (15) which ser v es as the inter mediate v elocity after enforcing the immersed-boundar y constraints. Bey ond correcting the velocity field, t he immersed-boundar y procedure also pro vides a direct route to evaluating the hydrodynamic force e xerted on the immersed body . At each sub-iteration 𝑚 , the direct-f orce term 𝐅 𝑚 𝐾 represents the momentum source per unit time required to satisfy t he no-slip condition at marker 𝐾 . Summing t he equal-and-opposite reaction over all markers and accumulating the contributions across all 𝑀 sub-iterations yields the net hydrodynamic f orce at the advanced time lev el: 𝐅 ( 𝑡 + Δ 𝑡 ) =  𝐾  − 𝑀 −1  𝑚 =0 𝐅 𝑚 𝐾 Δ 𝑠 𝐾  , (16) CL Li, H Xu: Preprint submitted to Elsevier P age 5 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws where the negativ e sign follo ws from N ewton ’ s third law , and the outer summation ov er sub-iterations accounts for the cumulative momentum impar ted to the fluid t hrough successiv e velocity cor rections (Eqs. 10 – 14 ). This f orce- ev aluation strategy is consistent with the multi-direct-forcing formulation of [ 42 ]. Subsequentl y , a pressure-correction step is applied to 𝐮 ∗ to recov er a div erg ence-free velocity consistent wit h the incompressibility constraint of the NS equations. In a conv entional projection framew ork, this cor rection is obt ained b y solving a pressure P oisson equation—an elliptic problem that requires repeated solutions of a large sparse linear system and typically dominates the computational cost of incompressible CFD. In t he proposed architecture (Fig. 2 ), t his bottlenec k is a voided by introducing a ResNe t-based implicit pressure-correction module. Specifically , t he diver gence field ∇ ⋅ 𝐮 ∗ (toge ther wit h the inter mediate v elocity f eatures) is provided as input to the Con vResNe t blocks, which learn to predict the pressure-induced cor rection and directly output the cor rected velocity at the next time lev el: 𝐮 ( 𝑡 + Δ 𝑡 ) =  𝑐 𝑜𝑛𝑣 2  𝐮 ∗ ( 𝑡 ) , ∇ ⋅ 𝐮 ∗ ( 𝑡 ); 𝜃 2  , (17) where 𝜃 2 is the trainable parameters of the ConvResN et module  𝑐 𝑜𝑛𝑣 2 . As a result, t he model enforces the projection effect without e xplicitly solving the pressur e Poisson system, yielding an efficient end-to-end differentiable surrogate while retaining the essential incompressibility constraint. 3.3. Single-step model training Previous studies ha ve ac hiev ed long-horizon stable inf erence using bac kpropagation t hrough time (BPTT)-based training sc hemes[ 43 ]; ho we ver , such strategies are no t well suited to ph ysics-integrated models. As repor ted in [ 36 ], BPTT training can be unstable in the early stages; moreov er, because t he ph ysics module is substantially more computationally demanding than the neural netw ork component, BPTT training typically incurs considerable GPU memory o verhead. In contrast, t he proposed differentiable physics-integrated networ k is trained using a one-step supervised strategy , where the model is optimized to matc h the ground-trut h flow fields and hydrodynamic f orces at the next time lev el without temporal unrolling. Given the input state at time 𝑡 , the model predicts the velocity components and the resultant f orce at 𝑡 + Δ 𝑡 , denoted by  𝑢 𝑡 +Δ 𝑡 ,  𝑣 𝑡 +Δ 𝑡 , and  𝐅 𝑡 +Δ 𝑡 , respectiv ely . The ov erall training objective is defined as:  ( 𝜃 1 , 𝜃 2 ) = 𝜆 𝐮 (  𝑢 ( 𝜃 1 , 𝜃 2 ) +  𝑣 ( 𝜃 1 , 𝜃 2 )) + 𝜆 𝐹  𝐅 ( 𝜃 1 ) , (18) where 𝜆 𝐮 and 𝜆 𝐅 weight the contr ibutions of the velocity and force super vision ter ms, respectivel y . Here we simply set 𝜆 𝐮 = 𝜆 𝐅 = 1 for all experiments to balance the two super vision ter ms. Con ventional data-nor malization heur istics do not, in general, conform to the pr inciple of phy sical similar ity (i.e., the scaling implied by similarity law s), which may lead to inconsistent er ror metr ics across heterogeneous targe ts. To ensure a consistent metr ic across velocity and for ce predictions, we adopt a relativ e 𝐿 𝑝 er ror (with 𝑝 = 2 in t his work) as the basic loss functional. For any prediction–targe t pair (  𝑦, 𝑦 ) , the relative loss is defined as:  𝑝 (  𝐲 , 𝐲 ) =   𝐲 − 𝐲  𝑝  𝐲  𝑝 + 𝜀 , (19) where  ⋅  𝑝 denotes the 𝐿 𝑝 norm computed over all degrees of freedom (i.e., after flattening t he cor responding field), and 𝜀 is a small constant f or numerical stability (set to 10 −8 in t his paper) to a void division by zero. Based on this relativ e 𝐿 𝑝 er ror , t he t hree loss components are compactly wr itten as:  𝑢 ( 𝜃 1 , 𝜃 2 ) = 1     𝑏 ∈   𝑝   𝑢 𝑡 +Δ 𝑡 𝑏 , 𝑢 𝑡 +Δ 𝑡 𝑏  ,  𝑣 ( 𝜃 1 , 𝜃 2 ) = 1     𝑏 ∈   𝑝   𝑣 𝑡 +Δ 𝑡 𝑏 , 𝑣 𝑡 +Δ 𝑡 𝑏   𝐅 ( 𝜃 1 ) = 1     𝑏 ∈   𝑝   𝐅 𝑡 +Δ 𝑡 𝑏 , 𝐅 𝑡 +Δ 𝑡 𝑏  , , (20) where  denotes t he mini-batc h index set. Here, 𝐅 𝑡 +Δ 𝑡 = [ 𝐹 𝑡 +Δ 𝑡 𝑥 , 𝐹 𝑡 +Δ 𝑡 𝑦 ] ⊤ denotes the ground-truth force v ector and  𝐅 𝑡 +Δ 𝑡 is the cor responding model prediction. CL Li, H Xu: Preprint submitted to Elsevier P age 6 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Finall y , t he model parameters are obt ained by minimizing the total loss ov er the training set via stoc hastic gradient descent: ( 𝜃 ∗ 1 , 𝜃 ∗ 2 ) = ar g min 𝜃 1 ,𝜃 2   ( 𝜃 1 , 𝜃 2 )  . (21) Because the optimization is per f or med on single-step targ ets, the training does not require bac kpropagation through long rollouts, which subs tantially reduces memory consumption and impr ov es training robustness f or ph ysics- integrated architectures. 4. Result and anal ysis 4.1. Baseline models for comparison T o better assess and benc hmark t he advantages of the proposed model, we consider tw o widely used classes of learning-based approaches: (i) a purely data-dr iven model and (ii) a physical loss-constrained model with additional phy sical loss ter ms. Schematic diagrams of these architectures are provided in Appendix C . In addition, we include a purely numerical baseline sol ved on the same coarse gr id and compare its results against the lear ning-based models, ensuring a comprehensive and fair evaluation framew ork. The purel y data-dr iv en netw ork, as shown in Fig. C2 (a), is constr ucted by removing all physics-based computation steps from our proposed model. Specifically , the fluid update is implemented using a ConvR esNet module with the same architecture as that emplo yed in our model. For t he structural dynamics component, to preserve a consistent data flow , t he original IB method computation block is replaced by an equiv alent ConvR esNet module, while all other components remain unc hanged. The total number of trainable parameters is kept comparable to t hat of our model, and we adopt the same training strategy and loss functions to eliminate confounding variables. Overall, the purel y data-dr iv en baseline consists of three independent R esNet ne tworks; its dat a flo w and training protocol are strictly aligned with those of our model, thereby enabling a fair and controlled comparison of per f or mance. For the phy sical loss-constrained model, as shown in Fig. C2 (b), the networ k architecture is identical to that of t he purely data-dr iven baseline. The only distinction lies in its loss function, whic h comprises a physics-based loss term and a data-fitting loss ter m, defined as follo ws:  t ot al = 𝛼  phy + 𝛽  𝑟 , (22) where 𝛼 and 𝛽 are w eighting coefficients selected follo wing the approach in [ 25 ] to ensure consistent ph ysical regularization strength. The data loss  𝑟 is kept identical to that used in the other models to ensure a f air compar ison; the ph ysics- consistency ter m  phy is specified as f ollow s:  phy =  𝑚 +  div +  IB , (23) where  𝑚 and  div cor respond to the momentum and continuity residual penalties, respectivel y , and  IB accounts for the boundar y -condition consistency at the immersed inter f ace.  𝑚 =      𝜕 𝐮 𝜕 𝑡 + ( 𝐮 ⋅ ∇) 𝐮 − 𝜈 ∇ 2 𝐮     2 2  , 𝐱 ∈ Ω 𝑓 ;  div =   ∇ ⋅ 𝐮  2  , 𝐱 ∈ Ω 𝑓 ;  IB =   𝑝 ( 𝐮 Γ , 𝐮 bc )  , 𝐱 ∈ Γ . (24) Here,  ⋅  denotes a veraging ov er the training batch and the corresponding spatiotemporal sample points (or grid locations), and Γ denotes the boundar y (and/or IB constraint) locations. In practice,  𝑚 and  div are ev aluated using a scheme-consistent discrete approximation that matches the numer ical update employ ed in the solv er, ensuring that the phy sical loss aligns with the underlying PDE discretization. CL Li, H Xu: Preprint submitted to Elsevier P age 7 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws 4.2. Flow past a stationary cir cular cy linder T o ev aluate the long-hor izon rollout capability of t he proposed model, we construct the training set o ver the nondimensional time window 𝑡 ∗ ∈ [0 , 100] , where 𝑡 ∗ = 𝑡𝑈 ∞ ∕ 𝐷 , and adopt t he same one-step supervised training strategy as described in Section 3 . After training, the model is initialized at 𝑡 ∗ = 0 and iterativel y rolled out until 𝑡 ∗ = 200 . Accordingly , predictions ov er 𝑡 ∗ ∈ [0 , 100] correspond to rollouts within the training hor izon, whereas those o ver 𝑡 ∗ ∈ [100 , 200] assess the model’ s e xtrapolation performance be yond the trained time range. The training hyperparameters are summar ized in Appendix B.2 f or reproducibility . W e compare t he proposed method against o ther baseline models from both spatial and temporal perspectives, as shown in Fig. 3 and Fig. 4 . Specificall y , the v or ticity snapshots assess t he models’ capability to infer complex spatial flow structures, whereas the time histories of str uctural loads ev aluate their long-horizon rollout st ability and accuracy . The results indicate that t he proposed model consistently outperforms all baseline models in both aspects, ev en though it is trained under a one-step supervised strategy using only non-sequential f eatures—highlighting the adv antage of embedding phy sical constraints into t he model architecture. By benc hmarking agains t the high-resolution numerical solution (the fifth column in Fig. 3 ), we obser v e that simulations at coarse resolution suffer from pronounced numer ical dissipation, leading to substantial prediction er rors. Moreov er, the ex cessive numer ical viscosity significantly delay s v or te x shedding relative to the high-resolution ref erence, demonstrating the limit ations of coarse-gr id numer ical solv ers for accurate flow prediction. In contras t, the purely data-dr iv en model div erg es rapidl y after only a short rollout horizon, which differs from the findings reported in [ 36 ]. This discrepancy mainly ar ises because, in the absence of embedded ph ysical constraints, the model lear ns a purely dat a-driven one-step evolution map and is then directly e xposed to long-horizon autoregressive rollout under single-step super vision. Under this setting, ev en small local prediction er rors are recursivel y f ed back into subsequent inputs, causing progressive dr ift of t he lear ned dynamics and ev entual breakdo wn. More impor tantly , suc h blac k - box autoregressiv e models do not explicitly enf orce the gov er ning ev olution structure or t he immersed-boundary interaction, so t he rollout can easil y depar t from t he phy sically admissible solution manifold once noticeable er rors appear . For this reason, man y e xisting CFD surrogate models[ 44 , 45 ] rely on explicitl y temporal arc hitectures (e.g., RNNs or Transf or mers) toge ther wit h BPTT training schemes to partially suppress error accumulation, albeit at the cost of substantially increased data requirements and training time. After incor porating the ph ysics-based loss term, the long-hor izon prediction accuracy of flow e volution is markedly impro ved, consis tent with pr ior research[ 46 ]. N ev er theless, compar ison with the g round truth rev eals a persistent phase drif t after prolonged rollouts, which is ev en more evident in s tr uctural force prediction (t he third ro w in Fig. 4 ). This beha vior indicates that soft physics regularization alone is insufficient to pre vent cumulative dynamical de viation: although t he loss ter m encourages physical consistency dur ing training, it does not explicitl y constrain t he rollout trajectory to follo w the numer ical update process or strictly preserve t he fluid–solid coupling at the solid boundary . As a result, boundary-condition mismatch and phase er ror can still accumulate over time, ev entually degrading both wake prediction and f orce estimation. In contrast, the proposed model accurately captures both t he near-cy linder flow field and the wak e dynamics from the onset of v or tex shedding to its full y dev eloped regime. This advantage is not merel y empirical, but stems directl y from t he model design. First, the discrete PDE operators are embedded into the residual update, so the network no long er learns an unconstrained blac k -box step-to-step mapping; instead, its e volution is restr icted by t he known numer ical phy sics, which effectivel y reduces the admissible solution space and suppresses the emergence of spur ious dynamics during long rollouts. Second, t he multi-resolution mechanism impro ves the accuracy of the immersed boundar y forcing, enabling more f ait hful resolution of the fluid–solid interaction and better enforcement of the boundary effect near the cy linder sur f ace. This is par ticularly import ant in the present problem, where long-ter m wak e de v elopment and f orce response are highly sensitiv e to t he accuracy of near-w all flow ev olution. Consequently , the proposed model show s substantially weaker error accumulation than both t he purel y dat a-driven model and the phy sics-loss-constrained model. A similar trend is obser v ed f or str uctural load prediction, as shown in the first ro w in Fig. 4 . Alt hough the earl y-stage dynamics f all within the training rollout range, the flow field in this regime remains highly transitional and t he vortex structures are still f or ming, making the f orce response especiall y sensitive to small phase and boundary er rors. As a result, both the purel y data-driven model and the ph ysics-loss-constrained model already accumulate larg e er rors at this stage. Benefiting from the explicit physical update pat h and t he more accurate treatment of immersed-boundary interaction, the proposed model maintains both amplitude and phase consis tency ov er long hor izons, producing structural load predictions that ag ree closely with the reference solution. CL Li, H Xu: Preprint submitted to Elsevier P age 8 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws T o facilitate a clearer compar ison of the long-hor izon rollout per f or mance across different models, t he relativ e er ror at each time step is defined as follo ws: 𝜖 𝑡 = 1 3   2   𝑢 𝑡 , 𝑢 𝑡  +  2   𝑣 𝑡 , 𝑣 𝑡  +  2   𝐅 𝑡 , 𝐅 𝑡   , (25) where 𝑡 denotes the time inde x of the rollout,  𝑢 𝑡 and  𝑣 𝑡 are t he predicted horizontal and v er tical velocity components at time 𝑡 , and 𝑢 𝑡 and 𝑣 𝑡 are t he cor responding ref erence solutions.  𝐅 𝑡 and 𝐅 𝑡 denote t he predicted and ref erence str uctural response v ectors, respectivel y . The operat or  2 ( ⋅ , ⋅ ) is the relativ e 𝐿 2 er ror defined in Eq. ( 19 ). The v elocity error is ev aluated separatel y f or the two components and then av eraged, whereas t he str uctural force er ror is computed using the full response vector as a whole. As shown in Fig. 5 , the purel y numerical solv er (black) exhibits a relative er ror that is sev eral orders of magnitude higher from the very beginning, pr imarily due to numerical viscosity and discretization er rors. This obser v ation furt her confir ms that, at coarse resolution, relying solely on low -fidelity numerical solv ers f or accelerated computation is no longer effective. The purel y data-dr iv en model (red) accumulates errors rapidly at the earl y stage of t he rollout and remains at a relativel y high lev el thereafter . The physics loss-constrained model (yello w) also show s rapid initial er ror g ro wth; the oscillations mainly stem from the phase mismatc h between the predictions and the ref erence solution. In contrast, the proposed model (blue) achiev es substantially smaller er rors than the af orementioned baselines throughout the rollout and effectiv ely suppresses error accumulation o ver long-horizon inf erence. 4.3. Flow past a rotationall y oscillating circular cylinder T o r igorously assess the performance of the proposed model and to mitigate potential bias arising from a single test scenario, we further evaluate t he model in the ro tationally oscillating cy linder configuration. The case setup follo ws that of the benchmark in Section 2 , ensur ing consistency wit h the over all problem f ormulation. T o better highlight the capability of the proposed approac h, we consider an angular -velocity amplitude of 𝜔 = 4 r ad∕s , cor responding to a non-lock -in regime. This regime is par ticularl y c hallenging because the wake dynamics are no long er gov erned by a single dominant shedding frequency; instead, t he y cont ain both intrinsic v or tex-shedding components and rotation- induced disturbances, resulting in weak er per iodicity and more intr icate fluid–str ucture interaction. Consequently , this case provides a more str ingent test of whe ther a surrogate model can preser ve physicall y consis tent long-hor izon ev olution bey ond simple periodic extrapolation. The training hyperparameters for all models are kept identical to those in Section 4.2 to ensure a fair comparison. The ev aluation results are presented in Figs. 6 and 7 , and the conclusions are larg ely consistent with those repor ted in Section 4.2 . At coarse resolution, the purely numerical solv er still substantially o veres timates the str uctural loads and exhibits considerable er rors in predicting t he vortex-shedding frequency . As the rollout horizon increases, both t he purely dat a-driven model and t he physics-loss-cons trained model e xhibit clear drift of the solid region, accompanied by progressiv ely amplified phase de viations in the wake v or tices. These er rors are more cr itical in the present non- lock -in regime, where the flow and f orce responses are less regular and therefore more sensitive to small inaccuracies in boundar y treatment and temporal ev olution. Once such er rors emerg e, t he y readily distort the coupled w ake dev elopment and str uctural loading, leading to cumulative degradation in long-horizon predictions. In contrast, the proposed physics-integrated neural differentiable model explicitl y enf orces t he solid-boundary flow beha vior through t he IB method and computes t he structural load response along a physicall y consistent update pat h. As a result, its adv antage in this case is not merel y higher instantaneous accuracy , but a stronger ability to maint ain stable and ph ysicall y admissible rollout trajectories in a more nonlinear and weakly per iodic regime. This is par ticularl y import ant because the non-lock -in condition r ules out t he possibility that good long-horizon performance is achiev ed simply b y reproducing an approximatel y per iodic pattern. Instead, the results indicate that the proposed model captures the underlying flo w-ev olution mechanism more faithfull y . As evidenced by the error cur ves in Fig. 8 , our model significantly outper f orms the other baselines and exhibits almost no error accumulation under long-hor izon rollouts, demonstrating super ior stability in a genuinel y challenging prediction scenar io. 5. Discussion 5.1. Generalizability ov er unseen cases In the e valuations repor ted in Section 4, all models perform temporal e xtrapolation under seen operating conditions. T o fur ther assess generalization to unseen conditions a cr itical cr iterion for practical sur rogate model deployment we CL Li, H Xu: Preprint submitted to Elsevier P age 9 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws consider the rotationally oscillating cylinder as an additional tes tbed. Specifically , the cases wit h rotation rates 𝜔 = 2 and 𝜔 = 4 are used for training, while 𝜔 = 3 and 𝜔 = 5 are reserved f or testing, representing an interpolative and an extrapolativ e regime, respectivel y . The training dat a cov er the nondimensional time window 𝑡 ∗ ∈ [0 , 100] , whereas f or the test conditions the trained model is rolled out ov er 𝑡 ∗ ∈ [0 , 200] to ev aluate both in-hor izon and out-of-horizon performance. By e xamining the mean rollout er ror computed as the relative er ror a veraged ov er the whole inf erence hor izon the results in Fig. 9 show t hat, in the cross-case comparison, the proposed model achiev es t he best per f or mance under both test conditions. In the wit hin-model comparison, all methods perform better on the inter polativ e case t han on the e xtrapolativ e case, which is consistent with the general expectation that inter polation is inherently easier than extrapolation. This is because the extrapolativ e condition is associated with a larger departure from the training regime, making the wak e ev olution and force response more difficult to reproduce accurately ov er long rollouts. T aking the more challenging e xtrapolative case with 𝜔 = 5 as an example, we furt her analyze the model performance from both temporal and spatial perspectives. The results are shown in Figs. 10 and 11 . The purel y data-dr iv en model and the ph ysics-loss-constrained model e xhibit behaviors larg ely consistent with t he conclusions in Section 4, namel y rapid er ror accumulation and eventual diver gence. By contrast, the proposed model sho ws only a moderate degradation under e xtrapolation: as the rollout horizon increases, the er rors in both the predicted flo w field and the s tr uctural load response gradually accumulate, which is more clearl y obser ved in F ig. 11 . Ne vertheless, benefiting from the constraints imposed by the ph ysicall y structured architecture, the proposed model still delivers substantially more accurate rollouts than the other approaches, even under this more demanding unseen condition. 5.2. T raining cost and inference acceleration ev aluation As shown in Fig. 12 , we repor t t he wall-cloc k time required for a single lear ning step Δ 𝑡 (cor responding to a time span of 100 numer ical time steps) and compare it across models. The purely data-dr iven model achiev es the shor test runtime, since it in vol ves no e xplicit ph ysics computation and relies solely on f ast ne twor k inf erence, resulting in the low est comput ational cost. The physics-loss-cons trained model shares the same netw ork architecture as the purely data-dr iv en baseline and is theref ore not shown separately . Owing to t he coarse gr id resolution and the larg er time step size, the coarse numer ical simulation is also substantially f aster t han t he high-resolution solv er, although it suffers from significant accuracy limit ations. The proposed model is slo wer than the purely data-driven baseline because it retains explicit physics-based computations wit hin the inference process. How ever , compared with conv entional numerical sol vers, it per mits a much larg er effective time step—f ar beyond the CFL-limited step size req uired by the embedded physics solv er itself—and a voids the expensive pressure-projection cor rection. As a result, on the same hardw are platform (a single RTX 3080Ti GPU with 12 GB memor y), the proposed approach achie ves appro ximately a 200-f old inference speedup relativ e to the high-resolution numerical solv er , and about a 20-f old speedup relative to the coarse-g rid numerical solv er at the same resolution, while maintaining good accuracy and phy sical consistency . In ter ms of training cost, we fur ther compare the single-step supervised strategy adopted in this wor k with backpr opagation through time (BPTT) using temporal windows of length 3 and 5. All e xperiments are conducted on the same personal works tation equipped wit h a single RTX 3080Ti GPU (12 GB memor y), using the same training hyperparameters to ensure a fair comparison. As sho wn in Fig. 13 , t he single-step strategy and BPTT with a 3-step window achie ve broadly comparable low -er ror per f or mance, whereas BPTT with a 5-step windo w shows slo wer con ver gence and remains at a noticeably higher er ror lev el wit hin the considered training time. This ma y be attributed to the limited model capacity and the finite amount of training data, whic h str uggle to suppor t the more complex optimization landscape of longer temporal windo ws. The competitiv e per f or mance ac hiev ed with only single-step training is mainly enabled by the phy sically structured design of the proposed model, where the main numerical update path and key physical cons traints are embedded directly into t he f or w ard propagation, thereby reducing t he dependence on long unrolled supervision f or learning stable temporal e volution. Consequentl y , within t he present framew ork, the model can achie ve competitiv e per f or mance using only single-step training, without implying that BPTT is generall y unnecessar y f or neural sur rogate modeling. A t the same time, the single-step setting is markedl y more efficient with respect to w all-clock training time, because eac h training update only in vol ves one-step f orward prediction and gradient computation, whereas BPTT must unroll the model o ver a temporal windo w and bac kpropagate through the cor responding computational graph. Even for the relativel y short window s considered here, this leads to substantially slow er con ver gence in ter ms of w all-clock time. CL Li, H Xu: Preprint submitted to Elsevier P age 10 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws T able 1 GPU memory consumption under different training strategies. T raining strategy GPU memory usage (GB) Single-step training 2.7 BPTT ( 𝑇 = 3 ) 5.7 BPTT ( 𝑇 = 5 ) 8.8 A similar trend is observed in memor y consumption. Because BPTT must retain intermediate activ ations ov er multiple rollout steps for gradient propagation, its memory usag e increases significantly with the temporal window length, as summar ized in Table 1 . In practice, this rapidly g ro wing memory demand is a major factor limiting the use of longer rollout horizons in training, and is one reason why many existing studies[ 39 , 40 ] resor t to tr uncated temporal window s to balance optimization q uality against computational cost. By contrast, t he proposed single-step training strategy av oids long-horizon graph construction altoge ther, enabling all cases in this work to be trained within one hour on a single consumer-grade GPU. This efficiency is particularl y valuable in settings where training dat a must be generated from expensiv e numer ical simulations and multiple rounds of model dev elopment are required. 6. Conclusion In this work, we ha ve de veloped a physics-integrated neural differentiable sur rogate f or incompressible flow s with immersed boundar ies, t arg eting accurate, stable, and efficient long-hor izon prediction. This work is an extension of the research by Fan and W ang, which fur ther advances the field of neural sur rogates f or fluid–str ucture interaction. The primar y contribution of this work lies in str ucturall y integrating physical pr inciples into the model architecture rather than treating physics as an auxiliary regular ization term. This integration is a ke y distinction from existing ph ysics- loss-constrained sur rogates. Specifically , t he proposed model embeds the pressure-projection update logic into an end- to-end differentiable framew ork, ensuring that t he model architecture and dat a flow are consistent with incompressible flow computation and t hat eac h intermediate variable retains a clear phy sical meaning. Within this framew ork, we make t hree additional k ey contr ibutions: an e xplicit immersed-boundary forcing module is incor porated to enf orce solid-boundary constraints and enable physicall y consistent hydrodynamic f orce evaluation; a sub-iteration strategy decouples the embedded phy sics solv er’ s internal stability requirement from t he sur rogate ’ s time step to enable stable coarse-gr id rollouts with larg e time s teps, which ov ercomes a cr itical limitation of conv entional coarse-grid solv ers; and a learned correction module replaces t he e xpensive pressure-P oisson projection to impro ve computational efficiency while preser ving accuracy . Notabl y , t he entire model can be trained via single-step supervision, which av oids long- horizon backpropagation and reduces training cost and memory consumption as another significant practical extension. Extensive tests on canonical benchmarks, including flow past a stationar y cylinder and a rotationally oscillating cylinder , v alidate our contr ibutions. These tests show that the proposed model consistently outperf or ms purel y data-dr iv en sur rogates, ph ysics-loss-constrained sur rog ates, and coarse-gr id numer ical sol vers in spatial fidelity and temporal robustness. The model accurately preserves near-body vortical structures, w ake ev olution, and str uctural load responses, while suppressing long-horizon error accumulation and phase dr ift. A furt her cr itical contribution is the validation of the model’ s generalization to unseen operating conditions. Using interpolative and extrapolativ e rotation rates, the model achie ves the low est mean rollout er ror, ev en under the more challenging extrapolativ e regime, which demonstrates broader applicability bey ond seen training conditions. More broadly , this w ork establishes a practical differentiable modeling framewor k for immersed-boundary flows. It e xtends the insights of Fan and W ang by demonstrating that physicall y structured arc hitecture design, rather than loss-based regularization alone, pr ovides a more effectiv e route tow ard stable long-hor izon prediction and efficient sur rogate simulation. These results confirm t hat the proposed approach is more than an accurate sur rog ate model. It delivers a novel, practical differentiable modeling framewor k for immersed-boundary flows, whic h is our o verarc hing contribution and a significant extension of Fan and W ang wit h clear potential f or scientific machine lear ning and differentiable simulation. By combining long-horizon stability , phy sical inter pretability , and high computational efficiency in a single trainable architectur e, the framew ork is well suited for applications requiring repeated flo w ev aluations, such as design optimization and flo w control, as it addresses a ke y comput ational bottlenec k. This advantage is par ticularl y rele vant f or emerging data-dr iven and reinf orcement-lear ning-based control paradigms, where the model can ser v e as an accurate CL Li, H Xu: Preprint submitted to Elsevier P age 11 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws and efficient environment sur rogate. Owing to its fully differentiable formulation and modular design, t he framew ork also pro vides a flexible foundation for future extensions to more complex fluid–str ucture interaction systems. A uthor contribution Chenglin Li: Conceptualization, Methodology , Sof tw are, For mal analysis, Inv estigation, W r iting – or iginal draf t. Hang X u: Conceptualization, Resources, Supervision, W riting – revie w & editing. Jianting Chen: Draft re vision. Y anfei Zhang: Draft revision. Data a vailability The data that suppor t the findings of this study are av ailable from the cor responding author upon request. A ckno wledgment This w ork was suppor ted by the open fund project of the State Ke y Laboratory of Maritime Tec hnology and Safe ty (Grant No.W25CG000074). Ref erences [1] J. Rabault, M. Kuchta, A. Jensen, U. Réglade, N. Cerardi, Ar tificial neural networks trained t hrough deep reinforcement learning discov er control strategies f or active flo w control, Journal of Fluid Mechanics 865 (2019) 281–302. doi: https://doi.org/10.1017/jfm.2019.62 . [2] S. Costanzo, T . Sayadi, M. Fosas de Pando, P . Schmid, P . Frey , Parallel-in-time adjoint-based optimization – application to unsteady incompressible flow s, Journal of Computational Physics 471 (2022) 111664. doi: https://doi.org/10.1016/j.jcp.2022.111664 . [3] W . Jia, H. X u, Strategies for energy -efficient flow control leveraging deep reinf orcement lear ning, Engineering Applications of Artificial Intelligence 156 (2025) 111125. doi: https://doi.org/10.1016/j.engappai.2025.111125 . [4] W . Jia, H. Xu, Optimal parallelization strategies for active flow control in deep reinforcement learning-based computational fluid dynamics, Phy sics of Fluids 36 (2024) 043623. doi: 10.1063/5.0204237 . [5] C. S. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 479–517. doi: 10.1017/S0962492902000077 . [6] J. F. Thompson, Z. U . W arsi, C. W ayne Mastin, Boundar y-fitted coordinate systems f or numerical solution of partial differential equations—a revie w , Jour nal of Computational Phy sics 47 (1982) 1–108. doi: https://doi.org/10.1016/0021- 9991(82)90066- 3 . [7] A . Johnson, T. T ezduyar, Mesh update strategies in parallel finite element computations of flow problems with moving boundar ies and in- terfaces, Computer Methods in Applied Mechanics and Engineering 119 (1994) 73–94. doi: https://doi.org/10.1016/0045- 7825(94) 00077- 8 . [8] R. Mittal, G. Iaccar ino, Immersed boundary methods, Annual Re view of Fluid Mechanics 37 (2005) 239–261. doi: https://doi.org/10. 1146/annurev.fluid.37.061903.175743 . [9] J. Y ang, F . Stern, A simple and efficient direct forcing immersed boundar y framew ork for fluid–str ucture interactions, Jour nal of Comput ational Phy sics 231 (2012) 5029–5061. doi: https://doi.org/10.1016/j.jcp.2012.04.012 . [10] H. Zhu, Y . Y ang, Z. Du, J. Y u, Gpu accelerated vortex-induced vibration simulation using jax: Efficiency and accuracy strategies, Computers & Fluids 305 (2026) 106913. doi: https://doi.org/10.1016/j.compfluid.2025.106913 . [11] A . E. Giannenas, S. Laizet, A simple and scalable immersed boundar y method for high-fidelity simulations of fixed and moving objects on a cartesian mesh, Applied Mathematical Modelling 99 (2021) 606–627. doi: https://doi.org/10.1016/j.apm.2021.06.026 . [12] R. Vinuesa, S. L. Brunton, Enhancing computational fluid dynamics wit h machine lear ning, Nature Comput ational Science 2 (2022) 358–366. doi: 10.1038/s43588- 022- 00264- 7 . [13] T . Davydzenka, P . T ahmasebi, High-resolution fluid–par ticle interactions: a mac hine lear ning approach, Journal of Fluid Mechanics 938 (2022) A20. doi: 10.1017/jfm.2022.174 . [14] K. T aira, S. L. Br unton, S. T . M. Da wson, C. W . Ro w ley , T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordey ev , V . Theofilis, L. S. Ukeiley , Modal analysis of fluid flows: An ov erview , AIAA Journal 55 (2017) 4013–4041. doi: 10.2514/1.J056060 . [15] T . Murata, K. Fukami, K. Fukagata, Nonlinear mode decomposition with conv olutional neural networ ks for fluid dynamics, Jour nal of Fluid Mechanics 882 (2020) A13. doi: 10.1017/jfm.2019.822 . [16] A . Cremades, S. Hoy as, R. Deshpande, P . Quintero, M. Lellep, W . J. Lee, J. P . Monty , N. Hutchins, M. Linkmann, I. Mar usic, R. Vinuesa, Identifying regions of impor tance in wall-bounded turbulence through explainable deep learning, Nature Communications 15 (2024) 3864. doi: 10.1038/s41467- 024- 47954- 6 . [17] T . Li, M. Buzzicotti, L. Bif erale, F . Bonaccorso, S. Chen, M. W an, Multi-scale reconstr uction of turbulent rotating flow s wit h proper or thogonal decomposition and generative adversarial networ ks, Journal of Fluid Mechanics 971 (2023) A3. doi: 10.1017/jfm.2023.573 . [18] F . Hadizadeh, W . Mallik, R. K. Jaiman, A graph neural networ k sur rogate model f or multi-objective fluid-acoustic shape optimization, Computer Methods in Applied Mechanics and Engineer ing 441 (2025) 117921. doi: https://doi.org/10.1016/j.cma.2025.117921 . [19] S. Barwey , H. Kim, R. Maulik, Interpretable a-posterior i er ror indication for graph neural network sur rogate models, Computer Methods in Applied Mechanics and Engineer ing 433 (2025) 117509. doi: https://doi.org/10.1016/j.cma.2024.117509 . [20] S. L. Brunton, B. R. N oack, P . Koumoutsak os, Machine lear ning f or fluid mechanics, Annual Revie w of Fluid Mec hanics 52 (2020) 477–508. doi: https://doi.org/10.1146/annurev- fluid- 010719- 060214 . CL Li, H Xu: Preprint submitted to Elsevier P age 12 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws [21] W . Liang, G. A. T adesse, D. Ho, L. Fei-Fei, M. Zahar ia, C. Zhang, J. Zou, A dvances, challeng es and oppor tunities in creating data for trustworthy ai, Nature Machine Intelligence 4 (2022) 669–677. doi: 10.1038/s42256- 022- 00516- 1 . [22] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P . Perdikaris, S. W ang, L. Y ang, Physics-inf or med machine lear ning, Nature Re view s Phy sics 3 (2021) 422–440. doi: 10.1038/s42254- 021- 00314- 5 . [23] J. T ompson, K. Schlachter , P . Sprechmann, K. Perlin, Accelerating euler ian fluid simulation with conv olutional networks, in: Proceedings of the 34th International Conf erence on Machine Learning - V olume 70, ICML ’17, JMLR.org, 2017, p. 3424–3433. [24] D. Koc hko v , J. A. Smith, A. Aliev a, Q. W ang, M. P . Brenner, S. Hoy er, Machine learning–accelerated comput ational fluid dynamics, Proceedings of the National Academy of Sciences 118 (2021) e2101784118. doi: 10.1073/pnas.2101784118 . [25] X.- Y . Liu, M. Zhu, L. Lu, H. Sun, J.-X. W ang, Multi-resolution partial differential equations preserved learning framew ork f or spatiotemporal dynamics, Communications Physics 7 (2024) 31. doi: 10.1038/s42005- 024- 01521- z . [26] X. Jin, S. Cai, H. Li, G. E. Kar niadakis, Nsfnets (navier -stokes flow nets): Physics-inf or med neural networ ks for the incompressible navier - stokes equations, Jour nal of Comput ational Physics 426 (2021) 109951. doi: https://doi.org/10.1016/j.jcp.2020.109951 . [27] L. Zhu, X. Jiang, A. Lefauv e, R. R. Kerswell, P . Linden, N ew insights into experimental stratified flows obtained through physics-inf ormed neural networks, Jour nal of Fluid Mechanics 981 (2024) R1. doi: 10.1017/jfm.2024.49 . [28] E. Haghighat, M. Raissi, A. Moure, H. Gomez, R. Juanes, A physics-inf or med deep learning framew ork for inv ersion and surrogate modeling in solid mechanics, Computer Methods in Applied Mechanics and Engineering 379 (2021) 113741. doi: https://doi.org/10.1016/j. cma.2021.113741 . [29] S. Cai, Z. W ang, S. W ang, P . Perdikaris, G. E. Kar niadakis, Physics-inf or med neural netw orks for heat transf er problems, Jour nal of Heat Transf er 143 (2021) 060801. doi: 10.1115/1.4050542 . [30] A . G. Baydin, B. A. Pearlmutter , A. A. Radul, J. M. Siskind, Automatic differentiation in machine lear ning: a sur v ey , Jour nal of Machine Learning Researc h 18 (2018) 1–43. [31] A . Krishnapriyan, A. Gholami, S. Zhe, R. Kirby , M. W . Mahone y , Characterizing possible failure modes in physics-inf ormed neural networks, in: M. Ranzato, A. Bey gelzimer , Y . Dauphin, P . Liang, J. W . V aughan (Eds.), A dvances in Neural Information Processing Systems, volume 34, Curran Associates, Inc., 2021, pp. 26548–26560. [32] R. T . Q. Chen, Y . Rubanov a, J. Bettencour t, D. Duvenaud, Neural ordinary differential equations, in: Proceedings of t he 32nd Inter national Conf erence on Neural Inf ormation Processing Systems, NIPS’18, Cur ran Associates Inc., Red Hook, NY , USA, 2018, p. 6572–6583. [33] Z. Long, Y . Lu, X. Ma, B. Dong, PDE-net: Learning PDEs from data, in: J. Dy , A. Krause (Eds.), Proceedings of t he 35th International Conf erence on Machine Learning, volume 80 of Proceedings of Machine Learning Researc h , PMLR, 2018, pp. 3208–3216. [34] X. Fan, X. Liu, M. W ang, J.-X. W ang, Diff-flowfsi: A gpu-optimized differentiable cfd platf orm f or high-fidelity turbulence and fsi simulations, 2025. . [35] D. A. Bezgin, A. B. Buhendwa, N. A. Adams, Jax-fluids: A fully-differentiable high-order computational fluid dynamics solv er f or compressible two-phase flo ws, Computer Ph ysics Communications 282 (2023) 108527. doi: https://doi.org/10.1016/j.cpc.2022.108527 . [36] X. Fan, J.-X. W ang, Differentiable hybrid neural modeling for fluid-str ucture interaction, Journal of Computational Phy sics 496 (2024) 112584. doi: https://doi.org/10.1016/j.jcp.2023.112584 . [37] X. Fan, D. Akhare, J.-X. W ang, Neural differentiable modeling wit h diffusion-based super-resolution f or two-dimensional spatiotemporal turbulence, Computer Methods in Applied Mechanics and Engineering 433 (2025) 117478. doi: https://doi.org/10.1016/j.cma.2024. 117478 . [38] D. Akhare, T. Luo, J.-X. W ang, Ph ysics-integrated neural differentiable (pindiff) model f or composites manufacturing, Computer Methods in Applied Mechanics and Engineer ing 406 (2023) 115902. doi: https://doi.org/10.1016/j.cma.2023.115902 . [39] W . W ang, X. Chu, Optimised flo w control based on automatic differentiation in compressible turbulent channel flow s, Journal of Fluid Mechanics 1011 (2025) A1. doi: 10.1017/jfm.2025.304 . [40] J. Xu, V . Mako viychuk, Y . N arang, F. Ramos, W . Matusik, A. Garg, M. Mac klin, Accelerated policy learning with parallel differentiable simulation, in: International Conf erence on Lear ning Representations, 2021, pp. 1–26. [41] S. Kumar, C. Lopez, O. Probst, G. Francisco, D. Askari, Y . Y ang, Flo w past a rotationally oscillating cylinder, Jour nal of Fluid Mechanics 735 (2013) 307–346. doi: 10.1017/jfm.2013.469 . [42] Z. W ang, J. Fan, K. Luo, Combined multi-direct f orcing and immersed boundary method f or simulating flows with moving particles, International Journal of Multiphase Flow 34 (2008) 283–302. doi: https://doi.org/10.1016/j.ijmultiphaseflow.2007.10.004 . [43] P . W erbos, Bac kpropagation through time: what it does and how to do it, Proceedings of the IEEE 78 (1990) 1550–1560. doi: 10.1109/5. 58337 . [44] K. Hasega wa, K. Fukami, T. Murat a, K. Fukagata, Cnn-lstm based reduced order modeling of two-dimensional unsteady flo ws around a circular cylinder at different reynolds numbers, Fluid Dynamics Research 52 (2020) 065501. doi: 10.1088/1873- 7005/abb91d . [45] Z. Li, T. Liu, W . Peng, Z. Y uan, J. W ang, A transformer -based neural operator f or large-eddy simulation of turbulence, Phy sics of Fluids 36 (2024) 065167. doi: 10.1063/5.0210493 . [46] C. Shi, G. Zhang, T . Sun, Y . Sui, Long-term prediction of cylinder water entr y flow fields using a physics-inf or med spatiotemporal graph neural network, Physics of Fluids 38 (2026) 023327. doi: 10.1063/5.0312013 . [47] X. Cui, X. Y ao, Z. W ang, M. Liu, A hybrid wa velet-based adaptive immersed boundary finite-difference lattice boltzmann method f or two- dimensional fluid–structure interaction, Journal of Computational Phy sics 333 (2017) 24–48. URL: https://www.sciencedirect.com/ science/article/pii/S0021999116306726 . doi: https://doi.org/10.1016/j.jcp.2016.12.019 . [48] M. Uhlmann, An immersed boundar y method with direct forcing f or the simulation of particulate flo ws, Journal of Computational Phy sics 209 (2005) 448–476. URL: https://www.sciencedirect.com/science/article/pii/S0021999105001385 . doi: https: //doi.org/10.1016/j.jcp.2005.03.017 . [49] Y . W ang, C. Shu, C. T eo, J. W u, An immersed boundary-lattice boltzmann flux solv er and its applications to fluid–structure interaction problems, Jour nal of Fluids and Str uctures 54 (2015) 440–465. URL: https://www.sciencedirect.com/science/article/pii/ CL Li, H Xu: Preprint submitted to Elsevier P age 13 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws S0889974614002709 . doi: https://doi.org/10.1016/j.jfluidstructs.2014.12.003 . Appendix A Numerical v erification A non-unif or m orthogonal Car tesian mesh with 717 × 382 cells is emplo yed o ver Ω 𝑓 . Grid stretc hing is applied a w ay from the body , while the mesh is locall y refined around t he cy linder within a rectangular region to resolv e near -wall dynamics; the minimum spacing near the cylinder is Δ 𝑥 min = Δ 𝑦 min = 𝐷 ∕64 . The immersed boundar y is represented by 𝑁 mar k er = 196 uniformly distributed Lagrangian markers along the cylinder sur face, and no-slip boundary conditions are imposed through direct forcing. T able A1 Dimensionless co efficients obtained from the simulation of the flow around a stationary cylinder at 𝑅𝑒 = 100 , using 𝐷 = 1 and 𝑑 𝑡 = 0 . 005 .  𝐶 𝐷 𝐶 ′ 𝐿 𝑆 𝑡 Our 1.383 ±0 . 345 0.167 Cui et al.[ 47 ] 1.360 ±0 . 340 0.167 Uhlmann[ 48 ] 1.453 ±0 . 339 0.169 W ang et al.[ 49 ] 1.334 ±0 . 370 0.163 Appendix B T raining details of neural model B.1 Downsam pling of high-resolution ground-truth data The training dat aset is constr ucted from the high-fidelity simulation data through two successive steps. First, in the spatial dimension, all field variables are transferred from t he fine gr id ( 717 × 382 cells ov er Ω 𝑓 ) to the coarse sur rogate gr id ( 80 × 40 cells ov er Ω 𝑟 𝑓 , cor responding to an 8× reduction in resolution) via the downsampling procedure illustrated in Fig. B1 : f ace-centered velocity components are area-a verag ed ov er each coarse control v olume, and cell-centered pressures are av eraged over the cor responding fine-gr id subcells. Second, in the temporal dimension, the spatiall y downsampled flow fields and structural forces are sampled at e very 100 high-fidelity time steps, yielding a learning time step of Δ 𝑡 = 100 × 𝑑 𝑡 = 0 . 5 . This sampling interval directly determines Δ 𝑡 , which defines t he prediction hor izon of the sur rogate model (see Eq. ( 7 )). B.2 Hyperparameters of the trainable Con vResNet block Each Con vResNe t block is composed of five con volutional la yers with channel configuration [32 , 64 , 64 , 32 , 32] and a trainable 3 × 3 ker nel. A rectified linear unit (ReLU) is applied af ter each la yer ex cept f or the las t la yer , where no activation is used to preserve the linear ity of t he output fields. The initial lear ning rate is set to 10 −3 , and the Adam optimizer is adopted with a weight decay of 10 −5 . A ReduceLR OnPlateau lear ning-rate scheduler is used, configured with mode = min, factor = 0.5, patience = 2, cooldo wn = 0, and min_lr = 10 −6 . The model is trained f or 3000 epochs with a batch size of 128. In our implementation, the lear ning-rate scheduler is driven by the validation loss rat her than the training loss. This choice is motivated by the intended deployment mode of the proposed model: it is ultimately used for autoregressive rollouts, where prediction er rors accumulate ov er time and the effectiv e per formance is gov er ned by generalization under distribution shift induced by t he model itself. Therefore, using the validation loss as the control signal provides a more reliable indicator of rollout-ready model quality than the one-s tep training objective alone, and yields a more sensible criter ion f or adapting the lear ning rate dur ing optimization. Appendix C Arc hitecture of baseline models(Fig. C2 ) CL Li, H Xu: Preprint submitted to Elsevier P age 14 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Figure 1: Schematic of the computational setup and domains. The immersed-boundary (IB) solver is defined on the global fluid domain Ω 𝑓 (gra y box) with boundaries Γ inlet , Γ out let , Γ upper , and Γ lower . The surrogate model is evaluated on a truncated sub domain Ω 𝑟 𝑓 (blue b o x) bounded by Γ 𝑟 inlet , Γ 𝑟 out let , Γ 𝑟 upper , and Γ 𝑟 lower . A prescribed inflow 𝑢 = 𝑢 in is imp osed at Γ inlet , and free-slip conditions are enforced on Γ upper and Γ lower . Figure 2: Schematic of the p rop osed physics-integrated differentiable neural netw ork. T rainable neural modules (ConvRes- Net blocks) a re shown in green, while non-trainable physics op erators (the PDE op erator and the multi-direct-forcing IBM mo dule) are highlighted in yello w. Grey dashed b oxes enclose the inputs, outputs, and intermediate physical va riables. Data propagate along the black a rrows. The PDE-operator substep is iterated 𝑁 times, and up-/down-sampling b ridges the resolutions b etw een the physics and neural components. CL Li, H Xu: Preprint submitted to Elsevier P age 15 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Figure 3: Comparison of spatio-temporal vo rticity predictions by different neural mo dels fo r flow past a stationa ry circular cylinder at ( 𝑅𝑒 = 100 , 𝜔 = 0) . Figure 4: Comparison of structural fo rce predictions by different neural mo dels for flow past a stationa ry circular cylinder at ( 𝑅𝑒 = 100 , 𝜔 = 0) . The left column sho ws the time histo ries of the drag coefficient 𝐶 𝐷 , and the right column shows the lift co efficient 𝐶 𝐿 . The gray-shaded interval denotes the training window, whereas the white interval co rresp onds to out-of-distribution time extrap olation b ey ond the training horizon. The black solid line rep resents the ground truth, and the blue curve co rresp onds to the prediction of the proposed mo del. CL Li, H Xu: Preprint submitted to Elsevier P age 16 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Figure 5: Stepwise rollout relative error 𝜖 𝑡 fo r flow past a stationary circula r cylinder at 𝑅𝑒 = 100 and 𝜔 = 0 , defined in Eq. ( 25 ) . The metric is computed as the average of the relative 𝐿 2 erro rs of 𝑢 , 𝑣 , and the structural resp onse vector 𝐅 at each rollout step. The gra y-shaded region indicates the training window, and the white region indicates temp oral extrap olation b ey ond the training horizon. Figure 6: Comparison of spatio-temporal vo rticity p redictions by different neural mo dels for flo w past a rotationally oscillating circular cylinder at ( 𝑅𝑒 = 100 , 𝜔 = 4) . CL Li, H Xu: Preprint submitted to Elsevier P age 17 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Figure 7: Compa rison of structural force p redictions by different neural mo dels for flow past a rotationally oscillating circular cylinder at ( 𝑅𝑒 = 100 , 𝜔 = 4) . The left column sho ws the time histo ries of the drag co efficient 𝐶 𝐷 , and the right column sho ws the lift co efficient 𝐶 𝐿 . The gray-shaded interval denotes the training window, whereas the white interval corresponds to out-of-distribution time extrap olation beyond the training horizon. The black solid line represents the ground truth, and the blue curve co rresp onds to the prediction of the proposed mo del. Figure 8: Stepwise rollout relative erro r 𝜖 𝑡 fo r flo w past a rotationally oscillating circular cylinder at 𝑅𝑒 = 100 and 𝜔 = 4 , defined in Eq. ( 25 ) . The metric is computed as the average of the relative 𝐿 2 erro rs of 𝑢 , 𝑣 , and the structural resp onse vecto r 𝐅 at each rollout step. The gray-shaded region indicates the training window, and the white region indicates temp oral extrap olation b ey ond the training horizon. CL Li, H Xu: Preprint submitted to Elsevier P age 18 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Figure 9: A verage relative erro r  𝜖 (averaged over the entire inference horizon) under unseen op erating conditions for the rotationally oscillating cylinder. Mo dels a re trained on 𝜔 = {2 , 4} and evaluated at 𝜔 = 3 (interp olation) and 𝜔 = 5 (extrap olation). Lo wer values indicate b etter long-horizon generalization. Figure 10: Compa rison of spatio-temporal vo rticity p redictions b y different neural mo dels for flo w past a rotationally oscillating circular cylinder at ( 𝑅𝑒 = 100 , 𝜔 = 5) . CL Li, H Xu: Preprint submitted to Elsevier P age 19 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Figure 11: Comparison of structural force p redictions by different neural mo dels for flo w past a rotationally oscillating circula r cylinder at ( 𝑅𝑒 = 100 , 𝜔 = 5) . The left column sho ws the time histo ries of the drag co efficient 𝐶 𝐷 , and the right column sho ws the lift co efficient 𝐶 𝐿 . The gray-shaded interval denotes the training windo w, whereas the white interval co rresp onds to out-of-distribution time extrap olation b ey ond the training horizon. The black solid line represents the ground truth, and the blue curve corresponds to the p rediction of the p rop osed mo del. Figure 12: W all-clo ck time p er learning step for different methods. The neural mo dels are evaluated on the truncated sub domain Ω 𝑟 𝑓 using an 80 × 40 grid, whereas the ground-truth CFD solution is computed on the full domain Ω 𝑓 using a 717 × 382 grid. CL Li, H Xu: Preprint submitted to Elsevier P age 20 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Figure 13: T raining efficiency compa rison betw een the single-step and BPTT-based multi-step training fo r flo w past a stationa ry circular cylinder at Reynolds number 𝑅𝑒 = 100 . The shaded bands represent ± 𝜎 . Figure B1: Illustration of transferring staggered-grid va riables from a fine mesh to an 𝑛 × 𝑛 coa rser mesh via downsampling (an 8× reduction in resolution is used in this study). Fo r each coarse control volume, the face-centered velo cit y comp onents on the fine grid a re a rea-averaged to obtain the corresponding coa rse-grid velo cities, while the cell-centered pressures within the fine-grid subcells are averaged to yield the coa rse-grid p ressure. Fine-grid velo city samples located in the interio r of the coa rse cell a re not retained. CL Li, H Xu: Preprint submitted to Elsevier P age 21 of 14 Physics-integrated neural differentiable modeling fo r IB flo ws Figure C2: Architectures of the baseline models used for compa rison. (a) Purely data-driven model trained with the data loss  𝑟 . (b) Physics loss constrained model, which augments (a) with the physics-consistency term  phy =  IB +  𝑚 +  div and optimizes the total loss  t ot al = 𝛼  phy + 𝛽  𝑟 . CL Li, H Xu: Preprint submitted to Elsevier P age 22 of 14