dragonfly_sim.core.time_integration =================================== .. py:module:: dragonfly_sim.core.time_integration .. autoapi-nested-parse:: Time integration: pseudo-transient continuation (PTC) for steady solves and BDF1/BDF2 for unsteady ones, on one residual. R*(U) = R(U) + inv_dt (a0 U + a1 U^n + a2 U^{n-1}) physical BDF term + (w_K / CFL) (U - U_old) pseudo-time term R is the model's steady residual (model.F, unchanged, which is what the adjoint differentiates). w_K = lambda_K / h_K (+ the model's time_scale_terms, e.g. viscous 2 nu / h^2) is the cell's inverse local time step at CFL 1, evaluated at U_old. Both terms are lumped DG0 mass terms and both live only in TimeTerm.augment(), so they never enter model.F. Modes (all one step loop, TimeIntegrator._newton_step): steady PTC physical term off (inv_dt = 0); one Newton step per pseudo step, U_old <- U before each, the CFL law (CFLController) adapts the step from the steady |R|. At the start of a pseudo step U = U_old, so the residual assembled there is the steady one. unsteady physical term on, pseudo term off (1/CFL = 0): Newton per physical step to the residual criterion; Delta t ramp and halving on failure (StepSizeController). dual time physical term on, and the pseudo term on inside every physical step, with the CFL law restarted each step: PTC on R* instead of R. Switching a term off sets its Constant to exactly zero, so the compiled forms stay valid across modes; a model with neither term enabled never builds them at all (the plain Newton solve is unchanged). .. !! processed by numpydoc !! Attributes ---------- .. autoapisummary:: dragonfly_sim.core.time_integration.BDF_COEFFS Classes ------- .. autoapisummary:: dragonfly_sim.core.time_integration.CFLController dragonfly_sim.core.time_integration.StepRecord dragonfly_sim.core.time_integration.StepSizeController dragonfly_sim.core.time_integration.TimeIntegrator dragonfly_sim.core.time_integration.TimeTerm Functions --------- .. autoapisummary:: dragonfly_sim.core.time_integration.bdf_coefficients Module Contents --------------- .. py:class:: CFLController Switched-evolution-relaxation CFL law with the safeguards the p = 0 RANS runs needed. Pure Python: update() takes the step's outcome and returns the decision, so the law is unit-tested without a solver. After a step from residual r_prev to r, with a converged linear solve (ksp_ok) and positivity step length theta: reject r non-finite or r > reject_factor r_prev: undo the step, CFL *= cut (stop when CFL < 1e-6 cfl0) limited theta < theta_cut: CFL *= cut, down to cfl_min_limited. Without it a limited step leaves |R| flat, the band below holds CFL, and the limiter goes on to zero every step (p=0 MUSCL from freestream froze at M 0.9 and 0.95; one leading-edge cell at p = 1e-4 blocked the step at CFL 1 but not at CFL 0.1, hence a floor below cfl_min). no Krylov an inexact direction never grows CFL progress r < (1 - 1e-3) r_prev: grow by max(r_prev/r, min_growth), capped at max_growth rise r > 1.05 r_prev: CFL *= r_prev/r, clipped to [rise_factor_min, rise_factor_max]; the defaults, 0.1 and 1/max_growth = 0.5, cut it by at least half band otherwise (|R| flat within -0.1%..+5%): grow by band_growth on a clean step (theta = 1), else hold. Holding kept the STW L3 wing at CFL 3.15 for 118 of 152 steps. CFL stays within [floor, cfl_max]; the floor is cfl_min_limited after a limited step and min(cfl_min, CFL) otherwise, so a CFL below cfl_min regrows by the usual rule instead of jumping back to cfl_min. exhausted: CFL fell below 1e-6 cfl0 after rejections, or the last max_frozen_steps steps all had theta = 0. A step with theta = 0 leaves the state unchanged, so once a cell sits at the positivity floor every later step repeats it; without this stop PTC spent its whole step budget frozen (OAT15A, Cadence grid level 1: ~950 identical steps at CFL cfl_min_limited). .. !! processed by numpydoc !! .. py:method:: __post_init__() .. py:method:: reset() .. py:method:: update(r_prev, r, ksp_ok=True, theta=1.0) Returns (accept, limited). Updates self.cfl. .. !! processed by numpydoc !! .. py:attribute:: band_growth :type: float :value: 1.2 .. py:attribute:: cfl :type: float .. py:attribute:: cfl0 :type: float :value: 5.0 .. py:attribute:: cfl_max :type: float :value: 1000000000000.0 .. py:attribute:: cfl_min :type: float :value: 1.0 .. py:attribute:: cfl_min_limited :type: float :value: 0.001 .. py:attribute:: cut :type: float :value: 0.1 .. py:property:: exhausted .. py:attribute:: frozen_steps :type: int :value: 0 .. py:attribute:: max_frozen_steps :type: int :value: 5 .. py:attribute:: max_growth :type: float :value: 2.0 .. py:attribute:: min_growth :type: float :value: 1.5 .. py:attribute:: reject_factor :type: float :value: 10.0 .. py:attribute:: rise_factor_max :type: float :value: 0.5 .. py:attribute:: rise_factor_min :type: float :value: 0.1 .. py:attribute:: theta_cut :type: float :value: 0.1 .. py:class:: StepRecord .. py:attribute:: accepted :type: bool .. py:attribute:: dt :type: float .. py:attribute:: ksp_its :type: int :value: 0 .. py:attribute:: newton_its :type: int :value: 0 .. py:attribute:: order :type: int :value: 0 .. py:attribute:: res_final :type: float .. py:attribute:: res_initial :type: float .. py:attribute:: retries :type: int :value: 0 .. py:class:: StepSizeController Physical time step: ramp from dt_start by `growth` per step up to dt, and halve on a failed step, up to max_retries times. .. !! processed by numpydoc !! .. py:method:: first(dt_prev) .. py:attribute:: dt :type: float .. py:attribute:: dt_start :type: float :value: None .. py:attribute:: growth :type: float :value: 1.5 .. py:attribute:: max_retries :type: int :value: 8 .. py:class:: TimeIntegrator(model, term) Drives a CompressibleEulerModel (or a model hosted on one) through steady PTC, unsteady BDF or dual-time steps. Built by the model (set_up_solver / enable_time_stepping); `term` is its TimeTerm. Time-dependent boundary conditions (core/farfield.py:TransverseFarfield) register through `boundary_updaters`: objects with prepare(order) before a physical step's Newton solve and commit() after it is accepted. .. !! processed by numpydoc !! .. py:method:: _copy(dst, src) .. py:method:: _dual_time_solve(ctrl, rtol, atol, max_it) PTC on R* inside one physical step; converged when |R*| < max(atol, rtol |R*_0|). Each pseudo step is one Newton step. .. !! processed by numpydoc !! .. py:method:: _newton_solve(rtol, atol, max_it) .. py:method:: _print(*args) .. py:method:: _unsteady_residual_norm(b) |R*| with the pseudo term off (at U = U_old it vanishes anyway). .. !! processed by numpydoc !! .. py:method:: _warm_start_cfl() The final CFL of the last converged steady solve if the state now in u_vec is that solve's end state (the warm start of an optimization step), else None. Collective. .. !! processed by numpydoc !! .. py:method:: advance(step_control, order=2, newton_rtol=1e-08, newton_atol=1e-10, newton_max_it=20, dual_time=None) One physical time step from U^n. The step size comes from step_control (StepSizeController); a failed Newton solve restores U^n and halves dt, up to step_control.max_retries times. `dual_time` (a CFLController) switches the inner solve to PTC on R*. Returns a StepRecord; record.accepted is False if every retry failed. .. !! processed by numpydoc !! .. py:method:: solve_steady(controller=None, max_steps=None, tol=None) Pseudo-transient continuation to the steady state from the current u_vec. Returns the model's physical residual Vec. With model.ptc_warm_start (default) and no explicit controller, a start from the end state of the previous converged steady solve continues at that solve's final CFL: the design (or alpha) has moved a little, so the flow is near its new steady state, and ramping the CFL up from cfl0 again would only repeat the early steps. A step that is too bold is rejected and the CFL cut, as always. .. !! processed by numpydoc !! .. py:method:: start(time=0.0, step_index=0, dt_prev=None) Begin (or restart) time stepping from the current u_vec, which is taken as U^n (and U^{n-1} unless restored from a checkpoint). .. !! processed by numpydoc !! .. py:attribute:: _steady_end :value: None .. py:attribute:: boundary_updaters :value: [] .. py:attribute:: dt_prev :value: None .. py:attribute:: model .. py:attribute:: step_index :value: 0 .. py:attribute:: term .. py:attribute:: time :value: 0.0 .. py:attribute:: verbose :value: True .. py:class:: TimeTerm(model, physical=False, pseudo=False) The physical and pseudo-time terms of R* as UFL forms over Constants and Functions owned here (see the module docstring). .. !! processed by numpydoc !! .. py:method:: augment(F) .. py:method:: local_inverse_time_scale(U) w_K = lambda_max / h_min plus the model's time_scale_terms. .. !! processed by numpydoc !! .. py:method:: physical_off() .. py:method:: set_bdf(order, dt, dt_prev=None) .. py:method:: set_cfl(cfl) .. py:attribute:: U_n .. py:attribute:: U_nm1 .. py:attribute:: U_old .. py:attribute:: a .. py:attribute:: inv_cfl .. py:attribute:: inv_dt .. py:attribute:: model .. py:function:: bdf_coefficients(order, ratio=None) (a0, a1, a2) for (a0 U^{n+1} + a1 U^n + a2 U^{n-1}) / dt_n. Variable-step BDF2 with ratio w = dt_n / dt_{n-1}: a0 = (1 + 2w)/(1 + w), a1 = -(1 + w), a2 = w^2/(1 + w), which is BDF_COEFFS[2] at w = 1. BDF1 ignores the ratio. .. !! processed by numpydoc !! .. py:data:: BDF_COEFFS