diff --git a/A_PLUS_APPROACHES_COMPARISON.md b/A_PLUS_APPROACHES_COMPARISON.md new file mode 100644 index 0000000..e6ac4f7 --- /dev/null +++ b/A_PLUS_APPROACHES_COMPARISON.md @@ -0,0 +1,82 @@ +# A+ (Plus Group) Implementation: Two Approaches Compared + +## Overview + +Both the `aplus` and `aplus3` branches implement the same fisheries "plus group" (A+) +accumulator — the oldest age class where survivors and newly-graduated individuals pool +together rather than continuing to age into new bins — but they take structurally different +approaches to scheduling and synchronizing A+ within the existing MPI cohort-parallel driver +(`main_cohort.cpp`). Both have been verified numerically equivalent: with A+ enabled they +produce identical checksums to each other, and with `-no-aplus` both exactly reproduce the +pre-A+ baseline, satisfying the conservation check `no-aplus checksum == (A+-on normal +checksum) + (A+-on A+ checksum)` to the last printed digit, on the real `skj_Fat.xml` config +(99/100 age groups, 480 months, 10 ranks). + +## `aplus`: graph-based scheduling + +A+ is folded directly into the same `SeapodymCohortDependencyAnalyzer` task graph that already +governs every ordinary cohort. Each calendar step of A+ becomes its own one-step task, appended +after the normal task IDs, with two explicit dependency edges: on A+'s own previous step, and on +the specific cohort step that "graduates" into it that month. Newborn cohorts get a third kind +of edge added automatically — a dependency on the matching A+ step, so a spawning cohort can +never be dispatched before A+'s density for that month is available. All of this is dispatched +by the same `TaskStepManager`/`TaskStepWorker` pair used for everything else; A+ tasks are just +another branch inside the existing `taskFunction`. No dedicated MPI rank, no extra communicator +split, no custom message protocol — correctness follows directly from the dependency graph's +existing "don't dispatch a task until every dependency is in the `completed` set" logic, a +mechanism that was already in use and already exercised long before A+ existed. + +## `aplus3`: dedicated worker with ping-pong handoff + +A+ is removed from the dependency graph entirely and instead runs as its own independent loop +on a dedicated MPI rank (`runAPlusWorker`). Cohorts that reach the oldest normal age class +("feeders") hand their graduating density directly to that rank via a blocking send, embedding +the payload in the message itself; A+ acknowledges once it has actually processed that calendar +step, and only after that ACK does the feeder tell the manager it's done — which is what +prevents anything downstream from running ahead of A+. Because feeders can complete out of +order, A+ carries a small reorder buffer that holds early arrivals and drains them strictly by +calendar step. This is a real custom protocol with its own state machine, and finding it correct +took real debugging effort — we found and fixed a genuine out-of-order feeding bug here that the +graph-based design was never exposed to, because it doesn't have a hand-rolled synchronization +protocol to get wrong in the first place. + +## Pros and cons + +| | `aplus` (graph) | `aplus3` (dedicated worker) | +|---|---|---| +| Minimum ranks | 2 | 3 | +| Extra machinery | None — reuses existing graph/dispatch | Dedicated worker loop, message tags, reorder buffer, ACK gating | +| Ordering guarantee | Declarative (graph edges) | Operational (blocking ACK + buffer) | +| A+ compute isolation | Shares dispatch queue with all cohorts | Own rank, never waits on farm backlog | +| Rank utilization | Every rank does real cohort work | A+'s rank is mostly idle (480 single-step advances vs. 50+ cohorts/rank elsewhere) | +| Debugging surface | Small — one well-tested mechanism | Larger — a protocol we had to design and validate ourselves | +| Verified correctness | Yes | Yes | + +`aplus`'s case is really about simplicity and efficiency: it adds no new synchronization +primitive, needs no extra rank, and every rank stays busy with real cohort work. Its downside is +that A+'s own advancement rides in the same queue as everything else, so if the farm is ever +fully saturated, A+ could in principle wait behind unrelated backlog — though the graph +guarantees it's never *wrong*, just potentially not the very next thing scheduled. + +`aplus3`'s case is about isolation: A+ gets guaranteed, uncontended throughput on its own rank, +decoupled from farm scheduling. The cost is a dedicated rank that, based on the timing captured +during testing, is doing meaningfully less work than any other rank (a single one-step-per-month +accumulator vs. dozens of full multi-step cohort trajectories), plus a bespoke protocol that's +more code to maintain and was the source of the one genuinely subtle bug found in this project +(the feeder-ordering gap). + +## Recommendation + +Converge on `aplus`. For the workload as it exists today — a single accumulator doing one step +of ordinary adult dynamics per calendar month — there's no evidence A+ is compute-bound enough +to need its own rank, and the graph-based design gets the same correctness guarantee with less +code, fewer ranks, and a synchronization mechanism that's already been proven out across the +rest of the scheduler rather than a new one built and debugged specifically for this feature. + +`aplus3` isn't wrong, and it's a legitimate design if A+'s per-step cost ever grows to actually +compete for farm time — a much heavier per-step computation, or multiple accumulator bins each +needing real throughput — but that's not the situation today, and paying for a dedicated rank +plus a bespoke protocol ahead of that need isn't buying anything. Worth keeping `aplus3` around +as a documented alternative rather than deleting it, since the two are now proven equivalent and +it remains a reasonable fallback if requirements change, but it shouldn't be the branch new work +builds on by default. diff --git a/example-configs/density_skj_fat/plot_workers_details.py b/example-configs/density_skj_fat/plot_workers_details.py index 660772f..6b27f67 100644 --- a/example-configs/density_skj_fat/plot_workers_details.py +++ b/example-configs/density_skj_fat/plot_workers_details.py @@ -61,7 +61,7 @@ def parse_logs(pattern): task_id = int(m.group(1)) task_ids.append(task_id) worker_ids.append(worker_id) - phases.append('coh init') + phases.append('init') t_starts.append(t_start) t_ends.append(ts) steps.append(-1) diff --git a/example-configs/density_skj_fat/rungenoa.sl b/example-configs/density_skj_fat/rungenoa.sl index 27b7428..6572b0d 100644 --- a/example-configs/density_skj_fat/rungenoa.sl +++ b/example-configs/density_skj_fat/rungenoa.sl @@ -20,4 +20,4 @@ module load gimkl/2020a rm log_*.txt # Run the MPI job -srun ../../bin/seapodym_cohort -s skj_fat.xml >& results/n${SLURM_NTASKS}.txt +srun ../../bin/seapodym_cohort_aplus -s skj_fat.xml >& results/n${SLURM_NTASKS}.txt diff --git a/example-configs/density_skj_fat/skj_fat.xml b/example-configs/density_skj_fat/skj_fat.xml index 7081781..d02bf59 100644 --- a/example-configs/density_skj_fat/skj_fat.xml +++ b/example-configs/density_skj_fat/skj_fat.xml @@ -244,7 +244,7 @@ - 1 2 47 + 1 2 97 @@ -255,7 +255,12 @@ 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 - 30 330 + 30 30 30 30 30 30 30 30 30 30 30 30 + 30 30 30 30 30 30 30 30 30 30 30 30 + 30 30 30 30 30 30 30 30 30 30 30 30 + 30 30 30 30 30 30 30 30 30 30 30 30 + 30 30 30 2970 + @@ -263,7 +268,12 @@ 48.33 51.19 53.86 56.36 58.70 60.88 62.92 64.83 66.61 68.27 69.83 71.28 72.64 73.91 75.10 76.21 77.25 78.22 79.12 79.97 80.76 81.50 82.19 82.83 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 - 85.00 87.96 + 85.00 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 + 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 + 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 + 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 + 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 + @@ -271,7 +281,13 @@ 2.27 2.73 3.21 3.72 4.23 4.76 5.30 5.83 6.36 6.89 7.40 7.91 8.41 8.89 9.36 9.81 10.25 10.66 11.07 11.45 11.82 12.17 12.51 12.83 13.5 13.5 13.5 13.5 13.50 13.50 13.50 13.50 13.50 13.50 13.50 13.50 - 14.0 15.56 + 14.0 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 + 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 + 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 + 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 + 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 + 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 + diff --git a/src/SeapodymCohort.cpp b/src/SeapodymCohort.cpp index 72f767d..7c50aff 100755 --- a/src/SeapodymCohort.cpp +++ b/src/SeapodymCohort.cpp @@ -3,6 +3,8 @@ #include "Date.h" #include "sys/stat.h" #include +#include +#include #include "DistDataCollector.h" #include "DataProvider.h" @@ -30,15 +32,15 @@ std::vector SeapodymCohort::GetCohortDensity() return vec; } -void SeapodymCohort::InitializeCohort(dvar_vector& x, DistDataCollector& dataCollector, const bool writeoutputfiles) +void SeapodymCohort::InitializeCohort(dvar_vector& x, DistDataCollector& dataCollector, int numTimeSteps, const bool writeoutputfiles) { - double t_all = MPI_Wtime(); - double t_rs = 0.0; - double t_reset = MPI_Wtime(); - //Reset model parameters: moved to worker init in main_cohort.cpp - //reset(x); - time_xreset += MPI_Wtime() - t_reset; +double t_all = MPI_Wtime(); +double t_rs = 0.0; +double t_reset = MPI_Wtime(); + //Reset model parameters: + reset(x); +time_xreset += MPI_Wtime() - t_reset; //----------------------------------------------// // ALLOCATE AND INITIALIZE COHORT DENSITY // @@ -88,6 +90,47 @@ double t_c = MPI_Wtime(); } t_copy_acc += MPI_Wtime() - t_c; } + + if (aPlusEnabled){ + // Plus group's density as of the end of t-1. SpawningBiomass_comp's + // own loop (a < param->sp_nb_cohorts[sp]) always included the + // oldest age class, A+ or not; nb_age_class here IS exactly that + // oldest (A+) index, since it excludes A+ from the "normal" range + // by construction. A+ publishes its density under this same + // chunk-id formula each step - see computeAPlusChunkId() and its + // use in main_cohort.cpp's taskFunction(). + int chunk_id = nb_age_class * numTimeSteps + (tstart_cohort - 1); + dataCollector.getAsync(chunk_id, data.data()); + +double t_f2 = MPI_Wtime(); + dataCollector.flush(); +t_flush_acc += MPI_Wtime() - t_f2; + + if (std::isnan(data[0])) { + // NaN (DistDataCollector::BAD_VALUE) here means chunk_id was + // never written before this cohort tried to read it - a + // missing/out-of-order dependency, not a numerical blow-up. + // Note: BAD_VALUE == NaN can never be detected with == + // (NaN != NaN by definition), hence std::isnan() here. + std::cerr << "ERROR: cohort " << cohort_id + << " (tstart_cohort=" << tstart_cohort << ") read BAD_VALUE (NaN) " + << "from A+ chunk " << chunk_id << " (t=" << (tstart_cohort - 1) << ")\n"; + MPI_Abort(MPI_COMM_WORLD, 1); + } + +double t_c2 = MPI_Wtime(); + int index = 0; + for (int i = map.imin1; i <= map.imax1; i++){ + const int jmin1 = map.jinf1[i]; + const int jmax1 = map.jsup1[i]; + for (int j = jmin1 ; j <= jmax1; j++){ + mat.dvarDensity(0,nb_age_class).elem_value(i,j) = data[index]; + index++; + } + } +t_copy_acc += MPI_Wtime() - t_c2; + } + dataCollector.endEpoch(); // should be as late as possible double t_total = MPI_Wtime() - t_block; @@ -119,6 +162,50 @@ double t_sp = MPI_Wtime(); time_spawning += MPI_Wtime() - t_sp; } + FinishInitialize(writeoutputfiles); + +time_init_cohort_restart += t_rs; +time_init_cohort_spawning += MPI_Wtime() - t_all - t_rs; +} + +void SeapodymCohort::InitializeAPlus(dvar_vector& x, const std::vector& mergedDensity, bool seedFromFile) +{ + //Reset model parameters: + reset(x); + + //----------------------------------------------// + // ALLOCATE AND INITIALIZE COHORT DENSITY // + //----------------------------------------------// + dvarCohortDensity.allocate(map.imin1, map.imax1, map.jinf1, map.jsup1); + if (seedFromFile){ + // t=0: no upstream task dependency; the plus group's opening balance + // is just the initial condition for its age bin (age_start == + // nb_age_class), read from file exactly like any other initial cohort. + RestoreDistributions(mat.nb_age_built); + dvarCohortDensity = mat.init_density_species(0,age_start); + } else { + // t>0: mergedDensity already holds (previous A+ pool) + (individuals + // that just graduated into the top age class this step), flattened + // in the same map.imin1..imax1/jinf1..jsup1 order GetCohortDensity() + // uses. This is the density stepForward() will apply this step's + // mortality/movement/feeding-habitat dynamics to - the AgePlus merge + // happens here, before dynamics, not instead of them. + int index = 0; + for (int i = map.imin1; i <= map.imax1; i++){ + const int jmin1 = map.jinf1[i]; + const int jmax1 = map.jsup1[i]; + for (int j = jmin1 ; j <= jmax1; j++){ + dvarCohortDensity.elem_value(i,j) = mergedDensity[index]; + index++; + } + } + } + + FinishInitialize(false); +} + +void SeapodymCohort::FinishInitialize(bool writeoutputfiles) +{ if (writeoutputfiles){ if (!param->gcalc()) ConsoleOutput(0,0); @@ -126,19 +213,19 @@ time_spawning += MPI_Wtime() - t_sp; //----------------------------------------------// // PRECOMPUTE VARIABLES USED IN COHORT MODELING // - //----------------------------------------------// + //----------------------------------------------// //precompute thermal habitat parameters for (int sp=0; sp < nb_species; sp++){ func.Vars_at_age_precomp(*param,sp); func.mortality_range_age_comp(*param,mat,sp); - + //precompute seasonal switch function if (param->seasonal_migrations[sp]){ func.Seasonal_switch_year_precomp(*param,mat,map, value(param->dvarsSpawning_season_peak[sp]), value(param->dvarsSpawning_season_start[sp]),sp); } - } + } getDate(jday, t_count);// In order to get qtr if (param->type_oxy==1) @@ -149,9 +236,6 @@ time_spawning += MPI_Wtime() - t_sp; age = age_start; past_month = month; past_qtr = qtr; - -time_init_cohort_restart += t_rs; -time_init_cohort_spawning += MPI_Wtime() - t_all - t_rs; } void SeapodymCohort::stepForward(bool writeoutputfiles) diff --git a/src/SeapodymCohort.h b/src/SeapodymCohort.h index ba98c0d..ef9aa3c 100755 --- a/src/SeapodymCohort.h +++ b/src/SeapodymCohort.h @@ -16,39 +16,66 @@ class SeapodymCohort : public SeapodymCoupled { public: SeapodymCohort(){/*DoesNothing*/}; - SeapodymCohort(const char* parfile, int cohortId) : SeapodymCoupled(parfile) { - cohort_id = cohortId; - nb_age_class = param->sp_nb_cohorts[0]; - - // Get starting age_class and start time from cohort_id - if (cohort_id >= nb_age_class){ - age_start = 0; - tstart_cohort = cohort_id-nb_age_class+1; - }else{ - age_start = nb_age_class - cohort_id - 1; - tstart_cohort = 0; - } + // aPlusOn selects whether the A+ (plus group) bin is carved out of the + // normal diagonal cohort scheme (default, matches main_cohort.cpp's + // default) or folded back in as an ordinary aging cohort, reproducing + // pre-A+ behavior. Exposed so main_cohort.cpp's -no-aplus flag can + // request the old behavior for regression comparison. + // nb_age_class excludes the A+ (plus group) bin when aPlusEnabled: the + // diagonal cohort-task scheme only ages "normal" cohorts through + // nb_age_class steps; the A+ bin (age index sp_nb_cohorts[0]-1) is + // handled separately as its own chain of accumulator tasks (see + // main_cohort.cpp / SeapodymCohortDependencyAnalyzer). When disabled, + // nb_age_class reverts to sp_nb_cohorts[0], i.e. A+ is just the last + // ordinary aging cohort, as before this feature existed. The + // cohort_id/age_start/tstart_cohort bookkeeping this implies is identical + // to what restart() computes, so delegate to it instead of duplicating + // it here - every worker's initial cohortId=0 is throwaway anyway (it + // gets restart()-ed to the real task id before any real work happens). + SeapodymCohort(const char* parfile, int cohortId, bool aPlusOn = true) : SeapodymCoupled(parfile) { + aPlusEnabled = aPlusOn; + restart(cohortId); }; virtual ~SeapodymCohort() {/*DoNothing*/}; - double time_overhead; - //double run_cohort(dvar_vector x, const bool writeoutputfiles = false) { return OnRunCohort(x, writeoutputfiles); } - void init_cohort(dvar_vector x, DistDataCollector& dataCollector, const bool writeoutputfiles = false) { return InitializeCohort(x, dataCollector, writeoutputfiles); } + double time_overhead; + //double run_cohort(dvar_vector x, const bool writeoutputfiles = false) { return OnRunCohort(x, writeoutputfiles); } + // numTimeSteps is only needed to locate A+'s published density chunk + // (aplus_chunk_id = nb_age_class*numTimeSteps + (tstart_cohort-1)) when + // this cohort is being initialized from spawning and aPlusEnabled - see + // InitializeCohort()'s spawning branch, and computeAPlusChunkId() below + // (same formula, expressed in terms of taskId/firstAPlusId instead of t). + void init_cohort(dvar_vector x, DistDataCollector& dataCollector, int numTimeSteps, const bool writeoutputfiles = false) { return InitializeCohort(x, dataCollector, numTimeSteps, writeoutputfiles); } void prerun_model(); void OnRunFirstStep(); double Checksum(); std::vector GetCohortDensity(); - int getChunkId(int step) { - int row = cohort_id - nb_age_class + 1 + step; + + // Chunk-id formula for a normal (non-A+) cohort task, factored out as a + // static so callers (e.g. main_cohort.cpp) can look up the chunk of *any* + // task/step pair, not just "this" cohort's own. + static int computeChunkId(int taskId, int step, int numAgeGroups) { + int row = taskId - numAgeGroups + 1 + step; int col = step; - if (cohort_idsp_nb_cohorts[0]; + nb_age_class = param->sp_nb_cohorts[0] - (aPlusEnabled ? 1 : 0); // see constructor's comment above // Get starting age_class and start time from cohort_id if (cohort_id >= nb_age_class){ age_start = 0; @@ -60,6 +87,36 @@ class SeapodymCohort : public SeapodymCoupled t_count = tstart_cohort+1; } + // Sets this object up to represent the A+ (plus group) accumulator task + // for absolute time index t (t=0 is the initial condition, before any + // time step has elapsed). age/age_start are pinned at nb_age_class (one + // past the last normal age, i.e. the A+ bin) rather than advancing: + // unlike restart(), each A+(t) is an independent, one-shot task/dispatch + // (see SeapodymCohortDependencyAnalyzer), not a persisting trajectory, so + // there is no "next age" to advance into. tstart_cohort/t_count follow + // the same convention restart() uses for a cohort "born" at time t, which + // is what makes stepForward()'s existing calendar-date and forcing-data + // lookups (both keyed off tstart_cohort/t_count) come out correct + // without any further changes to stepForward() itself. + void restartAPlus(int t){ + age_start = nb_age_class; + age = nb_age_class; + tstart_cohort = t; + t_count = t + 1; + } + + // Initializes the A+ task's density either from the actual initial + // condition (seedFromFile=true, used only for t=0, which has no upstream + // task dependency) or from mergedDensity - the previous A+ pool plus the + // individuals that just graduated into the top age class, already summed + // by the caller, flattened in the same map.imin1..imax1/jinf1..jsup1 + // order as GetCohortDensity(). Either way, stepForward() then runs the + // same adult dynamics (mortality, movement, feeding habitat) on it that + // any other cohort's step would get. + void init_cohort_aplus(dvar_vector x, const std::vector& mergedDensity, bool seedFromFile) { + return InitializeAPlus(x, mergedDensity, seedFromFile); + } + private: int dtau; int nbt_before_first_recruitment; @@ -78,6 +135,7 @@ class SeapodymCohort : public SeapodymCoupled int age_start; int tstart_cohort; int nb_age_class; + bool aPlusEnabled; DataProvider* dp_ = nullptr; @@ -98,7 +156,14 @@ class SeapodymCohort : public SeapodymCoupled int pop_built; - void InitializeCohort(dvar_vector& x, DistDataCollector& dataCollector, const bool writeoutputfiles = false); + void InitializeCohort(dvar_vector& x, DistDataCollector& dataCollector, int numTimeSteps, const bool writeoutputfiles = false); + void InitializeAPlus(dvar_vector& x, const std::vector& mergedDensity, bool seedFromFile); + // Setup shared by InitializeCohort() and InitializeAPlus(): precomputed + // per-age habitat/mortality parameters, calendar date, O2 climatology, + // and the age/past_month/past_qtr bookkeeping stepForward() relies on. + // Factored out so both initialization paths stay in lockstep instead of + // risking drift between two copies of the same ~15 lines. + void FinishInitialize(bool writeoutputfiles); public: void stepForward(const bool writeoutputfiles = false); diff --git a/src/main_cohort.cpp b/src/main_cohort.cpp index 42f2824..97cdc36 100755 --- a/src/main_cohort.cpp +++ b/src/main_cohort.cpp @@ -17,7 +17,6 @@ #include "Tags.h" #include #include -#include "admodel.h" double time_ic_comm = 0.0, time_ic_flush = 0.0, time_ic_copy = 0.0, time_spawning = 0.0, time_getdata = 0.0, time_xreset = 0.0, time_init_cohort_spawning = 0.0, time_init_cohort_restart = 0.0, time_io_forcing = 0.0; long n_ic = 0; // count of spawning-path inits, for per-init averages @@ -30,11 +29,11 @@ void buffers_set(long int &mv, long int &mc, long int &mg); double time_worker_init = 0.0, time_cohort_init = 0.0, time_calc = 0.0, time_mpi = 0.0, time_step = 0.0, time_overhead = 0.0; -SeapodymCohort xinit_prerun_wrapper(const char* parfile) { +SeapodymCohort xinit_prerun_wrapper(const char* parfile, bool useAPlus) { double tik = MPI_Wtime(); - SeapodymCohort cohort((char*)parfile, 0); + SeapodymCohort cohort((char*)parfile, 0, useAPlus); //initialize variables of optimization const int nvar = cohort.nvarcalc(); @@ -53,28 +52,98 @@ SeapodymCohort xinit_prerun_wrapper(const char* parfile) { return cohort; } +// Tell the manager task_id's given step is done. The third slot of the +// message used to be a separate "success" value that was always just a copy +// of task_id - collapsed here since it never carried any other information. +void notifyManagerDone(MPI_Comm comm, int task_id, int step) { + int output[3] = {task_id, step, task_id}; + MPI_Send(output, 3, MPI_INT, 0, END_TASK_TAG, comm); +} void taskFunction(int task_id, int stepBeg, int stepEnd, MPI_Comm comm, const std::shared_ptr& logger, DistDataCollector* dataCollector, SeapodymCohort* cohort, - const independent_variables& x){ + int firstAPlusId, int numAgeGroups, int numTimeSteps){ static double last_task_end = -1.0; // per-worker process, persists across calls double t_in = MPI_Wtime(); if (last_task_end >= 0.0) - time_idle += t_in - last_task_end; // <-- time spent in worker.run() waiting for dispatch - + time_idle += t_in - last_task_end; // <-- time spent in worker.run() waiting for dispatch + double tik = MPI_Wtime(); + //initialize variables of optimization - identical for both the A+ and + //normal-cohort paths below, so done once here rather than duplicated in + //each branch. + const int nvar = cohort->nvarcalc(); + independent_variables x(1, nvar); + adstring_array x_names(1,nvar); + cohort->xinit(x, x_names); + + if (task_id >= firstAPlusId) { + // A+ (plus group) accumulator task. Behaves like a normal cohort + // task from here on - initialize, then stepForward() runs the same + // mortality/movement/feeding-habitat dynamics any adult age gets - + // except its input density comes from merging two sources each step + // instead of a single spawning event, and its age is pinned rather + // than advancing (see SeapodymCohort::restartAPlus). stepBeg/stepEnd + // are always 0/1 for these tasks (see SeapodymCohortDependencyAnalyzer), + // so there is exactly one stepForward() call, not a loop. + int t = task_id - firstAPlusId; + logger->info("> A+ task id {} (t={})", task_id, t); + + cohort->restartAPlus(t); + + if (t == 0) { + // t=0: no upstream dependency; the plus group's opening balance + // is just the initial condition for its age bin, read from file + // like any other initial cohort. + cohort->init_cohort_aplus(x, std::vector(), /*seedFromFile=*/true); + } else { + int prevAPlusChunk = SeapodymCohort::computeAPlusChunkId(task_id - 1, firstAPlusId, numAgeGroups, numTimeSteps); + int graduatingChunk = SeapodymCohort::computeChunkId(t - 1, numAgeGroups - 1, numAgeGroups); + // Fetch the previous A+ pool directly into `merged`, then add the + // graduating cohort's density in place - avoids allocating a + // separate prevAPlus buffer just to sum it into a third one. + std::vector merged(dataCollector->getNumSize()); + std::vector graduating(merged.size()); + dataCollector->get(prevAPlusChunk, merged.data()); + dataCollector->get(graduatingChunk, graduating.data()); + for (std::size_t k = 0; k < merged.size(); ++k) + merged[k] += graduating[k]; + cohort->init_cohort_aplus(x, merged, /*seedFromFile=*/false); + } + + double tak_aplus = MPI_Wtime(); + time_cohort_init += tak_aplus - tik; + + // run this time step's real adult dynamics on the seeded/merged density + cohort->stepForward(false); + logger->info("End of A+ task id {} (t={}). Checksum = {}", task_id, t, cohort->Checksum()); + + std::vector out = cohort->GetCohortDensity(); + int myChunk = SeapodymCohort::computeAPlusChunkId(task_id, firstAPlusId, numAgeGroups, numTimeSteps); + + double t_p = MPI_Wtime(); + dataCollector->put(myChunk, out.data()); + time_mpi_put += MPI_Wtime() - t_p; + + notifyManagerDone(comm, task_id, stepBeg); + + time_calc += MPI_Wtime() - tak_aplus; + last_task_end = MPI_Wtime(); + logger->info("< A+ task id {} (t={})", task_id, t); + return; + } + logger->info("> task id {} for steps {} to {}", task_id, stepBeg, stepEnd); logger->info(" >> initialization of task id {}", task_id); - int cohort_id = task_id; cohort->restart(cohort_id); //initialize cohort either from restart or from spawning - cohort->init_cohort(x,*dataCollector); + cohort->init_cohort(x,*dataCollector,numTimeSteps); logger->info(" << initialization of task id {}", task_id); // advance the cohort @@ -100,11 +169,9 @@ void taskFunction(int task_id, int stepBeg, int stepEnd, MPI_Comm comm, time_mpi_put += MPI_Wtime() - t_p; logger->info(" <<< send data for step {} of task id {}", step, task_id); - int success = task_id; - int output[3] = {task_id, step, success}; logger->info(" >>> notify manager after step {} of task id {}", step, task_id); double tik_mpi = MPI_Wtime(); - MPI_Send(output, 3, MPI_INT, 0, END_TASK_TAG, comm); + notifyManagerDone(comm, task_id, step); time_mpi += MPI_Wtime() - tik_mpi; logger->info(" <<< notify manager after step {} of task id {}", step, task_id); } @@ -136,10 +203,13 @@ int main(int argc, char** argv) { CmdLineArgParser cmdLine; cmdLine.set("-s", std::string("initparfile.xml"), "Input parameter file"); + cmdLine.set("-no-aplus", false, "Disable the A+ (plus group) accumulator and reproduce " + "pre-A+ behavior/checksum, for regression comparison."); // Parse the command line arguments bool success = cmdLine.parse(argc, argv); bool help = cmdLine.get("-help") || cmdLine.get("-h"); + bool useAPlus = !cmdLine.get("-no-aplus"); if (!success) { std::cerr << "Error parsing command line arguments." << std::endl; cmdLine.help(); @@ -164,8 +234,14 @@ int main(int argc, char** argv) { param.init_param(); param.read(parfile); - // Get number of time steps and number of cohorts from param - int numAgeGroups = param.sp_nb_cohorts[0]; + // Get number of time steps and number of cohorts from param. + // numAgeGroups excludes the A+ (plus group) bin: the diagonal cohort-task + // scheme ages "normal" cohorts through numAgeGroups steps, and the A+ bin + // (age index sp_nb_cohorts[0]-1) is modelled separately as its own chain + // of one-step accumulator tasks (see SeapodymCohortDependencyAnalyzer). + // With -no-aplus, numAgeGroups reverts to sp_nb_cohorts[0] and A+ is just + // the last ordinary aging cohort, reproducing pre-A+ behavior. + int numAgeGroups = param.sp_nb_cohorts[0] - (useAPlus ? 1 : 0); int Tr_step, nbt_spinup_tuna, jday_run, jday_spinup, numTimeSteps; Date::init_time_variables(param, Tr_step, nbt_spinup_tuna, jday_run, jday_spinup, numTimeSteps, 0,0); //int numTasks = numAgeGroups + numTimeSteps - 1; @@ -177,22 +253,29 @@ int main(int argc, char** argv) { //Set-up the size for the shared arrays for forcing data std::vector> nameSizePairs = param.getDpNameSizePairs(numTimeSteps, map.get_array_size()); - // set up the data collector - int numChunks = numAgeGroups * numTimeSteps; + // set up the data collector. When A+ is enabled, one extra chunk per time + // step is reserved for the A+ (plus group) accumulator series, appended + // after the normal (task, step) chunk range - see + // SeapodymCohort::computeAPlusChunkId(). With -no-aplus this is 0, and the + // buffer is exactly the pre-A+ size. + int numChunks = numAgeGroups * numTimeSteps + (useAPlus ? numTimeSteps : 0); int color = (workerId == 0) ? 0 : 1; MPI_Comm workerComm; MPI_Comm_split(MPI_COMM_WORLD, color, workerId, &workerComm); if (workerId == 0) { - printf("[%d] Amount of data to be sent from workers to manager numData = %d numAgeGroups = %d numTimeSteps = %d numChunks = %d\n", \ - workerId, numData, numAgeGroups, numTimeSteps, numChunks); + printf("[%d] Amount of data to be sent from workers to manager numData = %d numAgeGroups = %d numTimeSteps = %d numChunks = %d aPlus = %s\n", \ + workerId, numData, numAgeGroups, numTimeSteps, numChunks, useAPlus ? "on" : "off"); } DistDataCollector dataCollect(MPI_COMM_WORLD, numChunks, numData); - // analyze the cohort Id task dependencies - SeapodymCohortDependencyAnalyzer taskDeps(numAgeGroups, numTimeSteps, param.age_mature[0]); + // analyze the cohort Id task dependencies (aPlusCohort adds the A+ chain; + // with -no-aplus this is false, and no A+ task ids are ever generated, + // so main()'s A+ branch below simply never triggers) + SeapodymCohortDependencyAnalyzer taskDeps(numAgeGroups, numTimeSteps, param.age_mature[0], /*aPlusCohort=*/useAPlus); + int firstAPlusId = taskDeps.getFirstAPlusCohortId(); int numCohorts = taskDeps.getNumberOfCohorts(); std::map stepBegMap = taskDeps.getStepBegMap(); std::map stepEndMap = taskDeps.getStepEndMap(); @@ -214,9 +297,18 @@ int main(int argc, char** argv) { // Make sure the data are ready for the final checksum MPI_Barrier(MPI_COMM_WORLD); double* data = dataCollect.getCollectedDataPtr(); - // print check sum - double checksum = std::accumulate(data, data + numChunks * numData, 0.0); - printf("[%d] Checksum = %15.5lf time manager = %10.5f sec\n", workerId, checksum, time_manager); + // Split the checksum into the "normal" cohort chunk range and the A+ + // chunk range (see SeapodymCohort::computeAPlusChunkId): normal chunks + // come first, A+ chunks are appended after. Keeping them separate lets + // us confirm a future change to A+'s dynamics doesn't perturb normal + // cohort results, and vice versa, instead of relying on one aggregate + // number that could mask a regression in either half. + int numNormalChunks = numAgeGroups * numTimeSteps; + double checksumNormal = std::accumulate(data, data + numNormalChunks * numData, 0.0); + double checksumAPlus = std::accumulate(data + numNormalChunks * numData, data + numChunks * numData, 0.0); + double checksum = checksumNormal + checksumAPlus; + printf("[%d] Checksum = %15.5lf (normal = %15.5lf, A+ = %15.8le) time manager = %10.5f sec\n", + workerId, checksum, checksumNormal, checksumAPlus, time_manager); //free manager's singleton 'color' MPI_Comm_free(&workerComm); @@ -248,15 +340,7 @@ int main(int argc, char** argv) { { DataProvider dp(workerComm, nameSizePairs); - SeapodymCohort cohort= xinit_prerun_wrapper(parfile.c_str()); - - // Initialize optimization variables once; x is stable for all tasks - // (xinit fills x from fixed model parameters that don't change between tasks). - // x is captured by reference in taskFunc — lifetime matches the enclosing block. - const int nvar = cohort.nvarcalc(); - independent_variables x(1, nvar); - adstring_array x_names(1, nvar); - cohort.xinit(x, x_names); + SeapodymCohort cohort= xinit_prerun_wrapper(parfile.c_str(), useAPlus); cohort.setDataProvider(&dp); @@ -266,11 +350,6 @@ int main(int argc, char** argv) { //MPI_Win_fence(0, dp.win()); double t_shm = MPI_Wtime(); cohort.setShmForcing(); // every node-local worker reads its timestep slice - - // Reset model parameters (applies boundp() transforms from optimisation space to - // physical parameter space) — done once here rather than per-task in InitializeCohort. - cohort.reset(dvar_vector(x)); - MPI_Barrier(dp.getShmComm()); // per-node publish-sync: all slabs visible before any read time_io_forcing += MPI_Wtime()-t_shm; @@ -283,7 +362,7 @@ int main(int argc, char** argv) { logger, &dataCollect, &cohort, - std::cref(x)); // x lives in this block, outlives taskFunc and worker.run() + firstAPlusId, numAgeGroups, numTimeSteps); TaskStepWorker worker(MPI_COMM_WORLD, taskFunc, stepBegMap, stepEndMap);