From c16fe9f167a19f21132d739bc173f6c190c90758 Mon Sep 17 00:00:00 2001 From: Sherwood Richers Date: Wed, 13 May 2026 14:21:24 -0400 Subject: [PATCH 1/4] first attempt at adaptive time integrator --- Src/Base/AMReX_IntegratorBase.H | 65 +++++++++++++++++++++++ Src/Base/AMReX_RKIntegrator.H | 91 +++++++++++++++++++++++++++++++-- Src/Base/AMReX_TimeIntegrator.H | 12 ++++- 3 files changed, 163 insertions(+), 5 deletions(-) diff --git a/Src/Base/AMReX_IntegratorBase.H b/Src/Base/AMReX_IntegratorBase.H index a6ba598aade..c5e7d07db6c 100644 --- a/Src/Base/AMReX_IntegratorBase.H +++ b/Src/Base/AMReX_IntegratorBase.H @@ -82,6 +82,45 @@ struct IntegratorOps } } + static amrex::Real ScaledMaxNorm (T& error, T& S_old, T& S_new, + amrex::Real abs_tol, amrex::Real rel_tol) { + // Compute max over all particles and real attributes of |error| / scale, + // where scale = abs_tol + rel_tol * max(|S_old|, |S_new|) per component. + // Using max(|S_old|, |S_new|) keeps the scale nonzero when a component crosses zero. + using TParIter = amrex::ParIter; + using ParticleType = amrex::Particle; + + amrex::ReduceOps reduce_op; + amrex::ReduceData reduce_data(reduce_op); + + int lev = 0; + TParIter pt_err(error, lev); + TParIter pt_old(S_old, lev); + TParIter pt_new(S_new, lev); + + for (; pt_err.isValid(); ++pt_err, ++pt_old, ++pt_new) { + const int np = pt_err.numParticles(); + const ParticleType* p_err = &(pt_err.GetArrayOfStructs()[0]); + const ParticleType* p_old = &(pt_old.GetArrayOfStructs()[0]); + const ParticleType* p_new = &(pt_new.GetArrayOfStructs()[0]); + reduce_op.eval(np, reduce_data, + [=] AMREX_GPU_DEVICE (int i) -> amrex::GpuTuple { + amrex::Real val = 0.0; + for (int c = 0; c < T::NStructReal; ++c) { + amrex::Real scale = abs_tol + rel_tol * amrex::max(std::abs(p_old[i].rdata(c)), + std::abs(p_new[i].rdata(c))); + val = amrex::max(val, std::abs(p_err[i].rdata(c)) / scale); + } + return {val}; + }); + } + + auto rv = reduce_data.value(); + amrex::Real result = amrex::get<0>(rv); + amrex::ParallelDescriptor::ReduceRealMax(result); + return result; + } + }; #endif @@ -130,6 +169,17 @@ struct IntegratorOps } } + static amrex::Real ScaledMaxNorm (T& error, T& S_old, T& S_new, + amrex::Real abs_tol, amrex::Real rel_tol) { + amrex::Real result = 0.0; + for (int i = 0; i < static_cast(error.size()); ++i) + for (int c = 0; c < error[i].nComp(); ++c) { + amrex::Real scale = abs_tol + rel_tol * amrex::max(S_old[i].norminf(c), S_new[i].norminf(c)); + result = amrex::max(result, error[i].norminf(c) / scale); + } + return result; + } + }; template @@ -160,6 +210,16 @@ struct IntegratorOps amrex::MultiFab::Saxpy(Y, a, X, scomp, scomp, mf_ncomp, nGrow); } + static amrex::Real ScaledMaxNorm (T& error, T& S_old, T& S_new, + amrex::Real abs_tol, amrex::Real rel_tol) { + amrex::Real result = 0.0; + for (int c = 0; c < error.nComp(); ++c) { + amrex::Real scale = abs_tol + rel_tol * amrex::max(S_old.norminf(c), S_new.norminf(c)); + result = amrex::max(result, error.norminf(c) / scale); + } + return result; + } + }; template @@ -340,6 +400,11 @@ public: return time_step; } + amrex::Real get_previous_time_step () + { + return previous_time_step; + } + void set_time_step (amrex::Real dt) { time_step = dt; diff --git a/Src/Base/AMReX_RKIntegrator.H b/Src/Base/AMReX_RKIntegrator.H index 4f9aebd4466..a881f8cd8a2 100644 --- a/Src/Base/AMReX_RKIntegrator.H +++ b/Src/Base/AMReX_RKIntegrator.H @@ -5,6 +5,7 @@ #include #include #include +#include namespace amrex { @@ -14,6 +15,7 @@ enum struct ButcherTableauTypes { Trapezoid, SSPRK3, RK4, + DormandPrince, NumTypes }; @@ -41,9 +43,22 @@ private: // RK embedded method b vector amrex::Vector extended_weights; + // Order of the embedded error estimate (e.g. 4 for Dormand-Prince) + int error_order; + + // Safety factor for adaptive step-size control + amrex::Real safety_factor; + + // Min/max allowed step-size growth/shrink factors per step + amrex::Real min_factor; + amrex::Real max_factor; + // RK stage right-hand sides amrex::Vector > F_nodes; + // Scratch space for low-order solution in embedded error estimate + amrex::Vector > S_error_vec; + // Current (internal) state and time amrex::Vector > S_current; amrex::Real time_current; @@ -84,6 +99,28 @@ private: {0.0, 0.0, 1.0, 0.0}}; weights = {1./6., 1./3., 1./3., 1./6.}; break; + case ButcherTableauTypes::DormandPrince: + // Dormand-Prince RK4(5) embedded pair + nodes = {0.0, + 1.0/5.0, + 3.0/10.0, + 4.0/5.0, + 8.0/9.0, + 1.0, + 1.0}; + tableau = {{0.0}, + {1.0/5.0, 0.0}, + {3.0/40.0, 9.0/40.0, 0.0}, + {44.0/45.0, -56.0/15.0, 32.0/9.0, 0.0}, + {19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0, -212.0/729.0, 0.0}, + {9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/18656.0, 0.0}, + {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0, 0.0}}; + // 5th-order weights + weights = {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0, 0.0}; + // 4th-order weights for embedded error estimate + extended_weights = {5179.0/57600.0, 0.0, 7571.0/16695.0, 393.0/640.0, -92097.0/339200.0, 187.0/2100.0, 1.0/40.0}; + error_order = 4; + break; default: amrex::Error("Invalid RK Integrator tableau type"); break; @@ -104,6 +141,16 @@ private: // By default, define no extended weights extended_weights = {}; + // Safety factor for adaptive step-size control (default 0.9) + safety_factor = 0.9; + pp.query("safety_factor", safety_factor); + + // Min/max allowed step-size growth/shrink factors per step + min_factor = 0.1; + max_factor = 5.0; + pp.query("min_factor", min_factor); + pp.query("max_factor", max_factor); + if (tableau_type == ButcherTableauTypes::User) { // Read weights/nodes/butcher tableau @@ -169,6 +216,11 @@ private: IntegratorOps::CreateLike(S_current, S_data, true); IntegratorOps::Copy(*S_current[0], S_data); + // Allocate error scratch space for embedded methods + if (!extended_weights.empty()) { + IntegratorOps::CreateLike(S_error_vec, S_data); + } + // Set the initial time time_current = time; } @@ -231,13 +283,46 @@ public: IntegratorOps::Saxpy(S_new, dt * weights[i], *F_nodes[i]); } - BaseT::post_step_action(S_new, time + dt); + BaseT::previous_time_step = dt; // If we are working with an extended Butcher tableau, we can estimate the error here, // and then calculate an adaptive time step. + amrex::Real dt_next = dt; + if (!extended_weights.empty() && BaseT::use_adaptive_time_step) + { + // Compute low-order solution: S_error = S_old + dt * sum(extended_weights[i] * F[i]) + IntegratorOps::Copy(*S_error_vec[0], S_old); + for (int i = 0; i < number_nodes; ++i) + IntegratorOps::Saxpy(*S_error_vec[0], dt * extended_weights[i], *F_nodes[i]); - // Save last completed step size for time_interpolate - BaseT::previous_time_step = dt; + // S_error = S_low - S_high (pointwise error estimate) + IntegratorOps::Saxpy(*S_error_vec[0], -1.0, S_new); + + const amrex::Real scaled_err = IntegratorOps::ScaledMaxNorm(*S_error_vec[0], S_old, S_new, + BaseT::abs_tol, BaseT::rel_tol); + + // Step-size controller: exponent is 1/(error_order+1) per standard embedded RK theory + amrex::Real factor = (scaled_err > 0.0) ? safety_factor * std::pow(1.0 / scaled_err, 1.0 / (error_order + 1)) : max_factor; + factor = amrex::min(amrex::max(factor, min_factor), max_factor); + dt_next = dt * factor; + + if (scaled_err > 1.0) + { + // Step rejected: retry with a smaller dt without calling post_step_action + return advance(S_old, S_new, time, dt_next); + } + + // Step accepted. Set time_step to max before calling post_step_action so that + // any CFL-based set_time_step call inside it can be detected and used to cap + // the adaptive suggestion. + BaseT::time_step = std::numeric_limits::max(); + } + + BaseT::post_step_action(S_new, time + dt); + + // Cap dt_next with any CFL limit set by post_step_action, then keep time_step in sync. + dt_next = amrex::min(dt_next, BaseT::time_step); + BaseT::time_step = dt_next; // Return time step return dt; diff --git a/Src/Base/AMReX_TimeIntegrator.H b/Src/Base/AMReX_TimeIntegrator.H index 2e1ce28e4fd..11554470be3 100644 --- a/Src/Base/AMReX_TimeIntegrator.H +++ b/Src/Base/AMReX_TimeIntegrator.H @@ -172,6 +172,11 @@ public: return integrator_ptr->get_time_step(); } + amrex::Real get_previous_time_step () + { + return integrator_ptr->get_previous_time_step(); + } + void set_time_step (amrex::Real dt) { integrator_ptr->set_time_step(dt); @@ -237,8 +242,11 @@ public: // Call the time integrator advance integrator_ptr->advance(S_old, S_new, m_time, m_timestep); - // Update our time variable - m_time += m_timestep; + // Update our time variable using the actual step size taken + m_time += integrator_ptr->get_previous_time_step(); + + // Pick up any timestep changes from post_step_action or adaptive control + m_timestep = integrator_ptr->get_time_step(); } } From 03e77207f21199024419bafc2cd8fc4b94caf46f Mon Sep 17 00:00:00 2001 From: Sherwood Richers Date: Sat, 23 May 2026 07:11:04 -0400 Subject: [PATCH 2/4] modifying error norm function, allowing for possibility of custom error norm, ability to recall previous error in adaptive integrator --- Src/Base/AMReX_IntegratorBase.H | 41 +++++++++++++++++++++++++-------- Src/Base/AMReX_RKIntegrator.H | 10 +++++--- Src/Base/AMReX_TimeIntegrator.H | 10 ++++++++ 3 files changed, 49 insertions(+), 12 deletions(-) diff --git a/Src/Base/AMReX_IntegratorBase.H b/Src/Base/AMReX_IntegratorBase.H index c5e7d07db6c..910cc16ce29 100644 --- a/Src/Base/AMReX_IntegratorBase.H +++ b/Src/Base/AMReX_IntegratorBase.H @@ -90,8 +90,7 @@ struct IntegratorOps using TParIter = amrex::ParIter; using ParticleType = amrex::Particle; - amrex::ReduceOps reduce_op; - amrex::ReduceData reduce_data(reduce_op); + amrex::Real result = amrex::Real(0.0); int lev = 0; TParIter pt_err(error, lev); @@ -100,23 +99,23 @@ struct IntegratorOps for (; pt_err.isValid(); ++pt_err, ++pt_old, ++pt_new) { const int np = pt_err.numParticles(); + if (np == 0) continue; const ParticleType* p_err = &(pt_err.GetArrayOfStructs()[0]); const ParticleType* p_old = &(pt_old.GetArrayOfStructs()[0]); const ParticleType* p_new = &(pt_new.GetArrayOfStructs()[0]); - reduce_op.eval(np, reduce_data, - [=] AMREX_GPU_DEVICE (int i) -> amrex::GpuTuple { + amrex::Real tile_max = amrex::Reduce::Max(np, + [=] AMREX_GPU_DEVICE (int i) -> amrex::Real { amrex::Real val = 0.0; for (int c = 0; c < T::NStructReal; ++c) { amrex::Real scale = abs_tol + rel_tol * amrex::max(std::abs(p_old[i].rdata(c)), std::abs(p_new[i].rdata(c))); val = amrex::max(val, std::abs(p_err[i].rdata(c)) / scale); } - return {val}; - }); + return val; + }, + amrex::Real(0.0)); + result = amrex::max(result, tile_max); } - - auto rv = reduce_data.value(); - amrex::Real result = amrex::get<0>(rv); amrex::ParallelDescriptor::ReduceRealMax(result); return result; } @@ -286,6 +285,14 @@ protected: */ std::function post_fast_step_action; + /** + * \brief Optional user-supplied error norm. If set, the adaptive + * step-size controller calls this instead of IntegratorOps::ScaledMaxNorm. + * Arguments are (error_state, S_old, S_new, abs_tol, rel_tol) and the + * return value is the scaled error used to decide step acceptance. + */ + std::function custom_error_norm; + /** * \brief Flag to enable/disable adaptive time stepping in single rate * methods or at the slow time scale in multirate methods (bool) @@ -346,6 +353,12 @@ protected: */ amrex::Real fast_abs_tol = 1.0e-9; + /** + * \brief Most recent scaled error computed by the adaptive controller. + * Updated each time advance() runs the error norm; otherwise zero. + */ + amrex::Real last_scaled_error = 0.0; + public: IntegratorBase () { @@ -395,6 +408,11 @@ public: post_fast_step_action = A; } + void set_error_norm (std::function F) + { + custom_error_norm = F; + } + amrex::Real get_time_step () { return time_step; @@ -405,6 +423,11 @@ public: return previous_time_step; } + amrex::Real get_scaled_error () + { + return last_scaled_error; + } + void set_time_step (amrex::Real dt) { time_step = dt; diff --git a/Src/Base/AMReX_RKIntegrator.H b/Src/Base/AMReX_RKIntegrator.H index a881f8cd8a2..57d77744938 100644 --- a/Src/Base/AMReX_RKIntegrator.H +++ b/Src/Base/AMReX_RKIntegrator.H @@ -298,11 +298,15 @@ public: // S_error = S_low - S_high (pointwise error estimate) IntegratorOps::Saxpy(*S_error_vec[0], -1.0, S_new); - const amrex::Real scaled_err = IntegratorOps::ScaledMaxNorm(*S_error_vec[0], S_old, S_new, - BaseT::abs_tol, BaseT::rel_tol); + const amrex::Real scaled_err = BaseT::custom_error_norm + ? BaseT::custom_error_norm(*S_error_vec[0], S_old, S_new, + BaseT::abs_tol, BaseT::rel_tol) + : IntegratorOps::ScaledMaxNorm(*S_error_vec[0], S_old, S_new, + BaseT::abs_tol, BaseT::rel_tol); + BaseT::last_scaled_error = scaled_err; // Step-size controller: exponent is 1/(error_order+1) per standard embedded RK theory - amrex::Real factor = (scaled_err > 0.0) ? safety_factor * std::pow(1.0 / scaled_err, 1.0 / (error_order + 1)) : max_factor; + amrex::Real factor = safety_factor * std::pow(1.0 / scaled_err, 1.0 / (error_order + 1)); factor = amrex::min(amrex::max(factor, min_factor), max_factor); dt_next = dt * factor; diff --git a/Src/Base/AMReX_TimeIntegrator.H b/Src/Base/AMReX_TimeIntegrator.H index 11554470be3..542e978e159 100644 --- a/Src/Base/AMReX_TimeIntegrator.H +++ b/Src/Base/AMReX_TimeIntegrator.H @@ -167,6 +167,11 @@ public: integrator_ptr->set_post_fast_step_action(A); } + void set_error_norm (std::function F) + { + integrator_ptr->set_error_norm(F); + } + amrex::Real get_time_step () { return integrator_ptr->get_time_step(); @@ -177,6 +182,11 @@ public: return integrator_ptr->get_previous_time_step(); } + amrex::Real get_scaled_error () + { + return integrator_ptr->get_scaled_error(); + } + void set_time_step (amrex::Real dt) { integrator_ptr->set_time_step(dt); From 36f8694d30df21f83dd114e8cfd6a08bd43b5ae4 Mon Sep 17 00:00:00 2001 From: Sherwood Richers Date: Tue, 2 Jun 2026 15:55:54 -0400 Subject: [PATCH 3/4] remove untested code --- Src/Base/AMReX_IntegratorBase.H | 63 ++------------------------------- Src/Base/AMReX_RKIntegrator.H | 9 +++-- 2 files changed, 6 insertions(+), 66 deletions(-) diff --git a/Src/Base/AMReX_IntegratorBase.H b/Src/Base/AMReX_IntegratorBase.H index 910cc16ce29..1d9e7eaf03e 100644 --- a/Src/Base/AMReX_IntegratorBase.H +++ b/Src/Base/AMReX_IntegratorBase.H @@ -82,44 +82,6 @@ struct IntegratorOps } } - static amrex::Real ScaledMaxNorm (T& error, T& S_old, T& S_new, - amrex::Real abs_tol, amrex::Real rel_tol) { - // Compute max over all particles and real attributes of |error| / scale, - // where scale = abs_tol + rel_tol * max(|S_old|, |S_new|) per component. - // Using max(|S_old|, |S_new|) keeps the scale nonzero when a component crosses zero. - using TParIter = amrex::ParIter; - using ParticleType = amrex::Particle; - - amrex::Real result = amrex::Real(0.0); - - int lev = 0; - TParIter pt_err(error, lev); - TParIter pt_old(S_old, lev); - TParIter pt_new(S_new, lev); - - for (; pt_err.isValid(); ++pt_err, ++pt_old, ++pt_new) { - const int np = pt_err.numParticles(); - if (np == 0) continue; - const ParticleType* p_err = &(pt_err.GetArrayOfStructs()[0]); - const ParticleType* p_old = &(pt_old.GetArrayOfStructs()[0]); - const ParticleType* p_new = &(pt_new.GetArrayOfStructs()[0]); - amrex::Real tile_max = amrex::Reduce::Max(np, - [=] AMREX_GPU_DEVICE (int i) -> amrex::Real { - amrex::Real val = 0.0; - for (int c = 0; c < T::NStructReal; ++c) { - amrex::Real scale = abs_tol + rel_tol * amrex::max(std::abs(p_old[i].rdata(c)), - std::abs(p_new[i].rdata(c))); - val = amrex::max(val, std::abs(p_err[i].rdata(c)) / scale); - } - return val; - }, - amrex::Real(0.0)); - result = amrex::max(result, tile_max); - } - amrex::ParallelDescriptor::ReduceRealMax(result); - return result; - } - }; #endif @@ -168,17 +130,6 @@ struct IntegratorOps } } - static amrex::Real ScaledMaxNorm (T& error, T& S_old, T& S_new, - amrex::Real abs_tol, amrex::Real rel_tol) { - amrex::Real result = 0.0; - for (int i = 0; i < static_cast(error.size()); ++i) - for (int c = 0; c < error[i].nComp(); ++c) { - amrex::Real scale = abs_tol + rel_tol * amrex::max(S_old[i].norminf(c), S_new[i].norminf(c)); - result = amrex::max(result, error[i].norminf(c) / scale); - } - return result; - } - }; template @@ -209,16 +160,6 @@ struct IntegratorOps amrex::MultiFab::Saxpy(Y, a, X, scomp, scomp, mf_ncomp, nGrow); } - static amrex::Real ScaledMaxNorm (T& error, T& S_old, T& S_new, - amrex::Real abs_tol, amrex::Real rel_tol) { - amrex::Real result = 0.0; - for (int c = 0; c < error.nComp(); ++c) { - amrex::Real scale = abs_tol + rel_tol * amrex::max(S_old.norminf(c), S_new.norminf(c)); - result = amrex::max(result, error.norminf(c) / scale); - } - return result; - } - }; template @@ -286,8 +227,8 @@ protected: std::function post_fast_step_action; /** - * \brief Optional user-supplied error norm. If set, the adaptive - * step-size controller calls this instead of IntegratorOps::ScaledMaxNorm. + * \brief User-supplied error norm. This must be set via set_error_norm + * before using adaptive time stepping; no default is provided. * Arguments are (error_state, S_old, S_new, abs_tol, rel_tol) and the * return value is the scaled error used to decide step acceptance. */ diff --git a/Src/Base/AMReX_RKIntegrator.H b/Src/Base/AMReX_RKIntegrator.H index 57d77744938..b545273d2e1 100644 --- a/Src/Base/AMReX_RKIntegrator.H +++ b/Src/Base/AMReX_RKIntegrator.H @@ -298,11 +298,10 @@ public: // S_error = S_low - S_high (pointwise error estimate) IntegratorOps::Saxpy(*S_error_vec[0], -1.0, S_new); - const amrex::Real scaled_err = BaseT::custom_error_norm - ? BaseT::custom_error_norm(*S_error_vec[0], S_old, S_new, - BaseT::abs_tol, BaseT::rel_tol) - : IntegratorOps::ScaledMaxNorm(*S_error_vec[0], S_old, S_new, - BaseT::abs_tol, BaseT::rel_tol); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(static_cast(BaseT::custom_error_norm), + "Adaptive time stepping requires an error norm to be set via set_error_norm."); + const amrex::Real scaled_err = BaseT::custom_error_norm(*S_error_vec[0], S_old, S_new, + BaseT::abs_tol, BaseT::rel_tol); BaseT::last_scaled_error = scaled_err; // Step-size controller: exponent is 1/(error_order+1) per standard embedded RK theory From c4a2ddb18bd7f549a83187892651a219cccc77c5 Mon Sep 17 00:00:00 2001 From: Sherwood Richers Date: Mon, 24 Aug 2026 12:07:04 -0400 Subject: [PATCH 4/4] add SoA capability to IntegratorOps::Saxpy --- Src/Base/AMReX_IntegratorBase.H | 37 ++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/Src/Base/AMReX_IntegratorBase.H b/Src/Base/AMReX_IntegratorBase.H index 1d9e7eaf03e..10028da42d4 100644 --- a/Src/Base/AMReX_IntegratorBase.H +++ b/Src/Base/AMReX_IntegratorBase.H @@ -44,8 +44,8 @@ struct IntegratorOps static void Saxpy (T& Y, const amrex::Real a, T& X) { // Calculate Y += a * X using a particle-level saxpy function supplied by the particle container T - using TParIter = amrex::ParIter; - using ParticleType = amrex::Particle; + using TParIter = typename T::ParIterType; + using ParticleType = typename T::ParticleType; int lev = 0; TParIter pty(Y, lev); @@ -69,16 +69,33 @@ struct IntegratorOps const int npx = ptx.numParticles(); AMREX_ALWAYS_ASSERT(npy == npx); - ParticleType* psy = &(pty.GetArrayOfStructs()[0]); - ParticleType* psx = &(ptx.GetArrayOfStructs()[0]); - auto particle_apply_rhs = T::particle_apply_rhs; - amrex::ParallelFor ( npy, [=] AMREX_GPU_DEVICE (int i) { - ParticleType& py = psy[i]; - const ParticleType& px = psx[i]; - particle_apply_rhs(py, a, px); - }); + // Hand particle_apply_rhs an array-of-structs particle when the + // container has one and the supplied function accepts it. + if constexpr ( ! ParticleType::is_soa_particle && + IsCallable::value) { + ParticleType* psy = &(pty.GetArrayOfStructs()[0]); + ParticleType* psx = &(ptx.GetArrayOfStructs()[0]); + + amrex::ParallelFor ( npy, [=] AMREX_GPU_DEVICE (int i) { + ParticleType& py = psy[i]; + const ParticleType& px = psx[i]; + particle_apply_rhs(py, a, px); + }); + + // Otherwise pass the particle tile data and an index, which is the only form + // available once the real components live in struct-of-arrays storage + } else { + auto ytd = pty.GetParticleTile().getParticleTileData(); + auto xtd = ptx.GetParticleTile().getParticleTileData(); + + amrex::ParallelFor ( npy, [=] AMREX_GPU_DEVICE (int i) { + particle_apply_rhs(ytd, i, a, xtd); + }); + } } }