rods/solver-c/tests/test_solver.c

#include "solver.h"
#include "solver_internal.h"

#include <math.h>
#include <stdio.h>
#include <string.h>

static void fill_base_inputs(SolverInputs *in) {
    memset(in, 0, sizeof(SolverInputs));
    in->schema_version = 2;
    in->workflow = 0;
    in->pumping_speed = 5.0;
    /* SI equivalents of typical imperial base-case magnitudes */
    in->pump_depth = 1727.0 * 0.3048;
    in->tubing_anchor_location = 1361.3 * 0.3048;
    in->rod_friction_coefficient = 0.2;
    in->stuffing_box_friction = 100.0 * 4.4482216152605;
    in->pump_friction = 200.0 * 4.4482216152605;
    in->taper_factor = 1.0;
    in->trajectory_friction_multiplier = 1.0;
    in->fluid_density_kg_m3 = 1000.0;
    in->gravity = 9.80665;
    in->upstroke_damping = 0.05;
    in->downstroke_damping = 0.15;
    in->non_dim_damping = 1.5;
    in->molded_guide_mu_scale = 1.5;
    in->wheeled_guide_mu_scale = 0.1;
    in->other_guide_mu_scale = 2.0;
}

static int test_deterministic_shape(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    SolverOutputs out1;
    SolverOutputs out2;

    if (solver_run(&inputs, &out1) != 0 || solver_run(&inputs, &out2) != 0) {
        return 1;
    }
    if (out1.point_count != 200 || out2.point_count != 200) {
        return 2;
    }
    for (int i = 0; i < out1.point_count; i++) {
        if (fabs(out1.polished_load[i] - out2.polished_load[i]) > 1e-9) {
            return 3;
        }
    }
    return 0;
}

static int test_bounds(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    SolverOutputs outputs;
    if (solver_run(&inputs, &outputs) != 0) {
        return 1;
    }
    if (fabs(outputs.max_polished_load - outputs.min_polished_load) < 1e-3) {
        return 2;
    }
    if (!isfinite(outputs.max_downhole_load) || !isfinite(outputs.min_downhole_load)) {
        return 3;
    }
    return 0;
}

static int test_fea_deterministic_shape(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    SolverOutputs out1;
    SolverOutputs out2;

    if (solver_run_fea(&inputs, &out1) != 0 || solver_run_fea(&inputs, &out2) != 0) {
        return 1;
    }
    if (out1.point_count != 200 || out2.point_count != 200) {
        return 2;
    }
    for (int i = 0; i < out1.point_count; i++) {
        if (fabs(out1.polished_load[i] - out2.polished_load[i]) > 1e-9) {
            return 3;
        }
    }
    return 0;
}

/* CFL diagnostic must be finite and physically plausible */
static int test_cfl_clamp(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    SolverOutputs outputs;
    if (solver_run_fdm(&inputs, &outputs) != 0) {
        return 1;
    }
    if (!(outputs.max_cfl > 0.0 && outputs.max_cfl < 1e6 && outputs.wave_speed_ref_m_s > 0.0)) {
        return 2;
    }
    return 0;
}

/* Static equilibrium helper: mean polished load should be finite and within broad physical band */
static int test_static_load_band(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    SolverOutputs outputs;
    if (solver_run_fdm(&inputs, &outputs) != 0) {
        return 1;
    }
    double sum = 0.0;
    for (int i = 0; i < outputs.point_count; i++) {
        sum += outputs.polished_load[i];
    }
    const double mean = sum / (double)outputs.point_count;
    if (!isfinite(mean) || mean <= 0.0) {
        return 2;
    }
    return 0;
}

/* Card extrema must match reported peaks (sanity / regression) */
static int test_card_extrema_match(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    SolverOutputs outputs;
    if (solver_run_fdm(&inputs, &outputs) != 0) {
        return 1;
    }
    double maxp = outputs.polished_load[0];
    double minp = outputs.polished_load[0];
    for (int i = 1; i < outputs.point_count; i++) {
        if (outputs.polished_load[i] > maxp) maxp = outputs.polished_load[i];
        if (outputs.polished_load[i] < minp) minp = outputs.polished_load[i];
    }
    if (fabs(maxp - outputs.max_polished_load) > 1e-3 || fabs(minp - outputs.min_polished_load) > 1e-3) {
        return 2;
    }
    return 0;
}

/* Undamped run remains finite (zero Rayleigh-like stroke factors) */
static int test_undamped_finite(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    inputs.upstroke_damping = 0.0;
    inputs.downstroke_damping = 0.0;
    inputs.non_dim_damping = 0.0;
    SolverOutputs outputs;
    if (solver_run_fdm(&inputs, &outputs) != 0) {
        return 1;
    }
    for (int i = 0; i < outputs.point_count; i++) {
        if (!isfinite(outputs.polished_load[i])) {
            return 2;
        }
    }
    return 0;
}

/* Zero harmonic drive: loads stay bounded (zero-input stability) */
static int test_zero_input_bounded(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    inputs.pumping_speed = 0.0;
    SolverOutputs outputs;
    if (solver_run_fdm(&inputs, &outputs) != 0) {
        return 1;
    }
    for (int i = 0; i < outputs.point_count; i++) {
        if (!isfinite(outputs.polished_load[i]) || fabs(outputs.polished_load[i]) > 1e12) {
            return 2;
        }
    }
    return 0;
}

/* Gas fraction and chamber pressure invariants after valve/gas stepping */
static int test_diagnostics_invariants(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    SolverOutputs outputs;
    if (solver_run_fdm(&inputs, &outputs) != 0) {
        return 1;
    }
    for (int i = 0; i < outputs.point_count; i++) {
        const double g = outputs.gas_fraction[i];
        if (g < -1e-6 || g > 1.0 + 1e-6) {
            return 2;
        }
        if (!isfinite(outputs.chamber_pressure_pa[i])) {
            return 3;
        }
    }
    return 0;
}

/* FDM vs FEA peak polished load tolerance (regression gate) */
static int test_fdm_fea_peak_tolerance(void) {
    SolverInputs inputs;
    fill_base_inputs(&inputs);
    SolverOutputs fdm;
    SolverOutputs fea;
    if (solver_run_fdm(&inputs, &fdm) != 0 || solver_run_fea(&inputs, &fea) != 0) {
        return 1;
    }
    const double tol = 700000.0; /* explicit FDM vs Newmark FEA — tighten after unified BC */
    if (fabs(fdm.max_polished_load - fea.max_polished_load) > tol) {
        return 2;
    }
    if (fabs(fdm.min_polished_load - fea.min_polished_load) > tol) {
        return 3;
    }
    return 0;
}

int main(void) {
    int rc = test_deterministic_shape();
    if (rc != 0) {
        printf("test_deterministic_shape failed: %d\n", rc);
        return 1;
    }

    rc = test_bounds();
    if (rc != 0) {
        printf("test_bounds failed: %d\n", rc);
        return 1;
    }

    rc = test_fea_deterministic_shape();
    if (rc != 0) {
        printf("test_fea_deterministic_shape failed: %d\n", rc);
        return 1;
    }

    rc = test_cfl_clamp();
    if (rc != 0) {
        printf("test_cfl_clamp failed: %d\n", rc);
        return 1;
    }

    rc = test_static_load_band();
    if (rc != 0) {
        printf("test_static_load_band failed: %d\n", rc);
        return 1;
    }

    rc = test_card_extrema_match();
    if (rc != 0) {
        printf("test_card_extrema_match failed: %d\n", rc);
        return 1;
    }

    rc = test_undamped_finite();
    if (rc != 0) {
        printf("test_undamped_finite failed: %d\n", rc);
        return 1;
    }

    rc = test_zero_input_bounded();
    if (rc != 0) {
        printf("test_zero_input_bounded failed: %d\n", rc);
        return 1;
    }

    rc = test_fdm_fea_peak_tolerance();
    if (rc != 0) {
        printf("test_fdm_fea_peak_tolerance failed: %d\n", rc);
        return 1;
    }

    rc = test_diagnostics_invariants();
    if (rc != 0) {
        printf("test_diagnostics_invariants failed: %d\n", rc);
        return 1;
    }

    printf("solver-c tests passed\n");
    return 0;
}