connermajic/rods
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Initial commit: establish deterministic rod-string solver stack.

Browse Source 
Set up the C solver core, Node API orchestration, TS GUI workflow, and engineering documentation with cleaned repo hygiene for private Git hosting.

Made-with: Cursor
This commit is contained in:
Conner Majic
2026-04-16 21:59:42 -06:00
commit 725a72a773
83 changed files with 14687 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
solver-c/CMakeLists.txt Normal file
View File

@@ -0,0 +1,28 @@
cmake_minimum_required(VERSION 3.16)
project(rod_solver_c C)
set(CMAKE_C_STANDARD 99)
set(CMAKE_C_STANDARD_REQUIRED ON)
add_library(rod_solver STATIC
src/solver_common.c
src/json_stdin.c
src/trajectory.c
src/solver_diagnostic.c
src/solver.c
src/solver_fea.c
src/solver_fourier.c
)
target_include_directories(rod_solver PUBLIC include)
add_executable(solver_main src/main.c)
target_link_libraries(solver_main PRIVATE rod_solver m)
add_executable(solver_fea_main src/main_fea.c)
target_link_libraries(solver_fea_main PRIVATE rod_solver m)
add_executable(test_solver tests/test_solver.c)
target_link_libraries(test_solver PRIVATE rod_solver m)
enable_testing()
add_test(NAME solver_c_tests COMMAND test_solver)

131
solver-c/include/solver.h Normal file
View File

@@ -0,0 +1,131 @@
#ifndef SOLVER_H
#define SOLVER_H
#include <stddef.h>
#define SOLVER_MAX_POINTS 512
#define SOLVER_MAX_WARNINGS 16
#define SOLVER_WARNING_LEN 160
#define SOLVER_MAX_NODES 65
#define SOLVER_MAX_SURVEY 512
#define SOLVER_MAX_SURFACE 4096
#define SOLVER_MAX_FOURIER_HARMONICS 32
typedef struct {
int schema_version;
int workflow; /* 0 predictive, 1 diagnostic */
double pumping_speed;
double pump_depth;
double tubing_anchor_location;
double rod_friction_coefficient;
double stuffing_box_friction;
double pump_friction;
double taper_factor;
double trajectory_friction_multiplier;
double fluid_density_kg_m3;
double gravity;
double upstroke_damping;
double downstroke_damping;
double non_dim_damping;
double molded_guide_mu_scale;
double wheeled_guide_mu_scale;
double other_guide_mu_scale;
int has_variable_rod;
int rod_node_count;
double area_m2[SOLVER_MAX_NODES];
double modulus_pa[SOLVER_MAX_NODES];
double density_kg_m3[SOLVER_MAX_NODES];
double weight_n_per_m[SOLVER_MAX_NODES];
double mts_n[SOLVER_MAX_NODES];
double guide_weight_n_per_m[SOLVER_MAX_NODES];
int survey_station_count;
double survey_md_m[SOLVER_MAX_SURVEY];
double survey_inc_rad[SOLVER_MAX_SURVEY];
double survey_azi_rad[SOLVER_MAX_SURVEY];
int geometry_valid;
double node_curvature[SOLVER_MAX_NODES];
double node_inc_rad[SOLVER_MAX_NODES];
double node_azi_rad[SOLVER_MAX_NODES];
double node_side_load_n[SOLVER_MAX_NODES];
int surface_count;
double surface_position_m[SOLVER_MAX_SURFACE];
double surface_load_n[SOLVER_MAX_SURFACE];
int surface_has_time;
double surface_time_s[SOLVER_MAX_SURFACE];
double pump_diameter_m;
double pump_intake_pressure_pa;
double tubing_id_m;
double percent_upstroke_time;
double percent_downstroke_time;
int pump_fillage_option;
double percent_pump_fillage;
double sinker_bar_diameter_m;
double sinker_bar_length_m;
int enable_profiles;
int enable_diagnostics_detail;
int enable_fourier_baseline;
int fourier_harmonics;
} SolverInputs;
typedef struct {
int point_count;
double position[SOLVER_MAX_POINTS];
double polished_load[SOLVER_MAX_POINTS];
double downhole_load[SOLVER_MAX_POINTS];
double pump_position_m[SOLVER_MAX_POINTS];
double pump_velocity_m_s[SOLVER_MAX_POINTS];
double polished_stress_pa[SOLVER_MAX_POINTS];
double side_load_profile_n[SOLVER_MAX_POINTS];
int profile_node_count;
double profile_md_m[SOLVER_MAX_NODES];
double profile_curvature_1pm[SOLVER_MAX_NODES];
double profile_inclination_rad[SOLVER_MAX_NODES];
double profile_azimuth_rad[SOLVER_MAX_NODES];
double profile_side_load_n[SOLVER_MAX_NODES];
double profile_friction_n[SOLVER_MAX_NODES];
int valve_traveling_open[SOLVER_MAX_POINTS];
int valve_standing_open[SOLVER_MAX_POINTS];
double chamber_pressure_pa[SOLVER_MAX_POINTS];
double gas_fraction[SOLVER_MAX_POINTS];
int fourier_harmonics_used;
double fourier_polished_load[SOLVER_MAX_POINTS];
double fourier_downhole_load[SOLVER_MAX_POINTS];
double fourier_residual_rms_polished;
double fourier_residual_rms_downhole;
double max_polished_load;
double min_polished_load;
double max_downhole_load;
double min_downhole_load;
int gas_interference;
double max_cfl;
double wave_speed_ref_m_s;
char warnings[SOLVER_MAX_WARNINGS][SOLVER_WARNING_LEN];
int warning_count;
} SolverOutputs;
/* json_stdin.c */
int solver_read_json_stdin(char **buffer_out, size_t *length_out);
int solver_parse_json_inputs(const char *buffer, SolverInputs *inputs);
/* trajectory.c */
void solver_trajectory_preprocess(SolverInputs *inputs, int nx, double rod_length_m);
/* solver_diagnostic.c */
int solver_run_diagnostic_fdm(const SolverInputs *inputs, SolverOutputs *outputs);
int solver_run_fdm(const SolverInputs *inputs, SolverOutputs *outputs);
int solver_run_fea(const SolverInputs *inputs, SolverOutputs *outputs);
int solver_run(const SolverInputs *inputs, SolverOutputs *outputs);
int solver_compute_fourier_baseline(const SolverInputs *inputs, SolverOutputs *outputs);
#endif

19
solver-c/include/solver_internal.h Normal file
View File

@@ -0,0 +1,19 @@
#ifndef SOLVER_INTERNAL_H
#define SOLVER_INTERNAL_H
#include "solver.h"
double solver_clamp(double v, double lo, double hi);
double solver_signum(double v);
void solver_add_warning(SolverOutputs *outputs, const char *msg);
void solver_init_output_ranges(SolverOutputs *outputs);
void solver_update_output_ranges(SolverOutputs *outputs, double polished, double downhole);
double solver_input_or_default(double value, double fallback);
double solver_compute_side_load_node(const SolverInputs *inputs, double tension_n, int node_idx, double ds);
double solver_compute_friction_node(const SolverInputs *inputs, double side_load_n, double velocity_m_s, int node_idx);
void solver_fill_profiles(const SolverInputs *inputs, SolverOutputs *outputs, int node_count, double rod_length_m,
const double *side_load_nodes, const double *friction_nodes);
void solver_valve_state_step(const SolverInputs *inputs, SolverOutputs *outputs, int step_idx, double pump_position_m,
double pump_velocity_m_s, double downhole_load_n);
#endif

360
solver-c/src/json_stdin.c Normal file
View File

@@ -0,0 +1,360 @@
#include "solver.h"
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define JSON_READ_CHUNK 65536
#define JSON_MAX_SIZE (4 * 1024 * 1024)
int solver_read_json_stdin(char **buffer_out, size_t *length_out) {
size_t cap = JSON_READ_CHUNK;
size_t len = 0;
char *buf = (char *)malloc(cap);
if (!buf) {
return -1;
}
for (;;) {
size_t n = fread(buf + len, 1, JSON_READ_CHUNK, stdin);
len += n;
if (n < JSON_READ_CHUNK) {
break;
}
if (len + JSON_READ_CHUNK > JSON_MAX_SIZE) {
free(buf);
return -2;
}
if (len + JSON_READ_CHUNK > cap) {
cap *= 2;
char *nb = (char *)realloc(buf, cap);
if (!nb) {
free(buf);
return -3;
}
buf = nb;
}
}
buf[len] = '\0';
*buffer_out = buf;
*length_out = len;
return 0;
}
static void skip_ws(const char **p) {
while (**p && isspace((unsigned char)**p)) {
(*p)++;
}
}
static const char *find_matching_brace(const char *p) {
if (*p != '{') {
return NULL;
}
int depth = 0;
const char *s = p;
for (; *s; s++) {
if (*s == '{') depth++;
else if (*s == '}') {
depth--;
if (depth == 0) {
return s;
}
}
}
return NULL;
}
static int parse_string_literal(const char **p, char *out, size_t out_cap) {
skip_ws(p);
if (**p != '"') {
return -1;
}
(*p)++;
size_t i = 0;
while (**p && **p != '"' && i + 1 < out_cap) {
out[i++] = *(*p)++;
}
out[i] = '\0';
if (**p == '"') {
(*p)++;
}
return 0;
}
static int parse_number(const char **p, double *out) {
skip_ws(p);
char *end = NULL;
double v = strtod(*p, &end);
if (end == *p) {
return -1;
}
*p = end;
*out = v;
return 0;
}
static int parse_bool(const char **p, int *out) {
skip_ws(p);
if (strncmp(*p, "true", 4) == 0) {
*p += 4;
*out = 1;
return 0;
}
if (strncmp(*p, "false", 5) == 0) {
*p += 5;
*out = 0;
return 0;
}
return -1;
}
static int parse_double_array(const char **p, double *out, int max_n, int *count_out) {
skip_ws(p);
if (**p != '[') {
return -1;
}
(*p)++;
int n = 0;
for (;;) {
skip_ws(p);
if (**p == ']') {
(*p)++;
break;
}
if (n >= max_n) {
return -2;
}
if (parse_number(p, &out[n]) != 0) {
return -3;
}
n++;
skip_ws(p);
if (**p == ',') {
(*p)++;
continue;
}
if (**p == ']') {
(*p)++;
break;
}
return -4;
}
*count_out = n;
return 0;
}
static const char *object_find(const char *obj_start, const char *obj_end, const char *key) {
char pattern[96];
snprintf(pattern, sizeof(pattern), "\"%s\"", key);
const char *hit = obj_start;
while (hit < obj_end) {
hit = strstr(hit, pattern);
if (!hit || hit >= obj_end) {
return NULL;
}
if (hit > obj_start && (isalnum((unsigned char)hit[-1]) || hit[-1] == '_')) {
hit += 1;
continue;
}
const char *colon = strchr(hit + strlen(pattern), ':');
if (!colon || colon >= obj_end) {
return NULL;
}
colon++;
while (colon < obj_end && isspace((unsigned char)*colon)) {
colon++;
}
return colon;
}
return NULL;
}
static int parse_object_number(const char *obj_start, const char *obj_end, const char *key, double *out) {
const char *p = object_find(obj_start, obj_end, key);
if (!p) {
return -1;
}
const char *q = p;
if (parse_number(&q, out) != 0) {
return -2;
}
return 0;
}
static int parse_object_int(const char *obj_start, const char *obj_end, const char *key, int *out) {
double v = 0.0;
if (parse_object_number(obj_start, obj_end, key, &v) != 0) {
return -1;
}
*out = (int)v;
return 0;
}
static int parse_object_bool(const char *obj_start, const char *obj_end, const char *key, int *out) {
const char *p = object_find(obj_start, obj_end, key);
if (!p) {
return -1;
}
const char *q = p;
return parse_bool(&q, out);
}
static int parse_object_array(const char *obj_start, const char *obj_end, const char *key, double *arr, int max_n, int *count_out) {
const char *p = object_find(obj_start, obj_end, key);
if (!p) {
return -1;
}
const char *q = p;
return parse_double_array(&q, arr, max_n, count_out);
}
int solver_parse_json_inputs(const char *buffer, SolverInputs *inputs) {
if (!buffer || !inputs) {
return -1;
}
memset(inputs, 0, sizeof(SolverInputs));
inputs->schema_version = 2;
inputs->workflow = 0;
inputs->gravity = 9.80665;
inputs->fluid_density_kg_m3 = 1000.0;
inputs->taper_factor = 1.0;
inputs->trajectory_friction_multiplier = 1.0;
inputs->molded_guide_mu_scale = 1.0;
inputs->wheeled_guide_mu_scale = 1.0;
inputs->other_guide_mu_scale = 1.0;
inputs->percent_upstroke_time = 50.0;
inputs->percent_downstroke_time = 50.0;
inputs->enable_profiles = 0;
inputs->enable_diagnostics_detail = 0;
inputs->enable_fourier_baseline = 0;
inputs->fourier_harmonics = 8;
const char *wf = strstr(buffer, "\"workflow\"");
if (wf) {
const char *colon = strchr(wf, ':');
if (colon) {
colon++;
while (*colon && isspace((unsigned char)*colon)) colon++;
char tmp[32];
const char *q = colon;
if (parse_string_literal(&q, tmp, sizeof(tmp)) == 0) {
if (strcmp(tmp, "diagnostic") == 0) {
inputs->workflow = 1;
}
}
}
}
const char *model_key = strstr(buffer, "\"model\"");
if (!model_key) {
return -2;
}
const char *brace = strchr(model_key, '{');
if (!brace) {
return -3;
}
const char *model_end = find_matching_brace(brace);
if (!model_end) {
return -4;
}
/* required scalars with defaults if missing */
(void)parse_object_number(brace, model_end, "pumping_speed", &inputs->pumping_speed);
(void)parse_object_number(brace, model_end, "pump_depth", &inputs->pump_depth);
(void)parse_object_number(brace, model_end, "tubing_anchor_location", &inputs->tubing_anchor_location);
(void)parse_object_number(brace, model_end, "rod_friction_coefficient", &inputs->rod_friction_coefficient);
(void)parse_object_number(brace, model_end, "stuffing_box_friction", &inputs->stuffing_box_friction);
(void)parse_object_number(brace, model_end, "pump_friction", &inputs->pump_friction);
(void)parse_object_number(brace, model_end, "taper_factor", &inputs->taper_factor);
(void)parse_object_number(brace, model_end, "trajectory_friction_multiplier", &inputs->trajectory_friction_multiplier);
(void)parse_object_number(brace, model_end, "fluid_density_kg_m3", &inputs->fluid_density_kg_m3);
(void)parse_object_number(brace, model_end, "gravity", &inputs->gravity);
(void)parse_object_number(brace, model_end, "upstroke_damping", &inputs->upstroke_damping);
(void)parse_object_number(brace, model_end, "downstroke_damping", &inputs->downstroke_damping);
(void)parse_object_number(brace, model_end, "non_dim_damping", &inputs->non_dim_damping);
(void)parse_object_number(brace, model_end, "molded_guide_mu_scale", &inputs->molded_guide_mu_scale);
(void)parse_object_number(brace, model_end, "wheeled_guide_mu_scale", &inputs->wheeled_guide_mu_scale);
(void)parse_object_number(brace, model_end, "other_guide_mu_scale", &inputs->other_guide_mu_scale);
(void)parse_object_number(brace, model_end, "pump_diameter_m", &inputs->pump_diameter_m);
(void)parse_object_number(brace, model_end, "pump_intake_pressure_pa", &inputs->pump_intake_pressure_pa);
(void)parse_object_number(brace, model_end, "tubing_id_m", &inputs->tubing_id_m);
(void)parse_object_number(brace, model_end, "percent_upstroke_time", &inputs->percent_upstroke_time);
(void)parse_object_number(brace, model_end, "percent_downstroke_time", &inputs->percent_downstroke_time);
(void)parse_object_int(brace, model_end, "pump_fillage_option", &inputs->pump_fillage_option);
(void)parse_object_number(brace, model_end, "percent_pump_fillage", &inputs->percent_pump_fillage);
(void)parse_object_number(brace, model_end, "sinker_bar_diameter_m", &inputs->sinker_bar_diameter_m);
(void)parse_object_number(brace, model_end, "sinker_bar_length_m", &inputs->sinker_bar_length_m);
int hr = 0;
if (parse_object_bool(brace, model_end, "has_variable_rod", &hr) == 0) {
inputs->has_variable_rod = hr;
}
(void)parse_object_int(brace, model_end, "rod_node_count", &inputs->rod_node_count);
int na = 0, ne = 0, nd = 0;
if (parse_object_array(brace, model_end, "area_m2", inputs->area_m2, SOLVER_MAX_NODES, &na) == 0) {
(void)na;
}
if (parse_object_array(brace, model_end, "modulus_pa", inputs->modulus_pa, SOLVER_MAX_NODES, &ne) == 0) {
(void)ne;
}
if (parse_object_array(brace, model_end, "density_kg_m3", inputs->density_kg_m3, SOLVER_MAX_NODES, &nd) == 0) {
(void)nd;
}
int nw = 0, nm = 0, ng = 0;
(void)parse_object_array(brace, model_end, "weight_n_per_m", inputs->weight_n_per_m, SOLVER_MAX_NODES, &nw);
(void)parse_object_array(brace, model_end, "mts_n", inputs->mts_n, SOLVER_MAX_NODES, &nm);
(void)parse_object_array(brace, model_end, "guide_weight_n_per_m", inputs->guide_weight_n_per_m, SOLVER_MAX_NODES, &ng);
if (na > 0 && ne > 0 && nd > 0 && na == ne && na == nd) {
inputs->has_variable_rod = 1;
inputs->rod_node_count = na;
}
int ns = 0;
if (parse_object_array(brace, model_end, "survey_md_m", inputs->survey_md_m, SOLVER_MAX_SURVEY, &ns) == 0 && ns > 0) {
int ni = 0, nz = 0;
if (parse_object_array(brace, model_end, "survey_inc_rad", inputs->survey_inc_rad, SOLVER_MAX_SURVEY, &ni) == 0 &&
parse_object_array(brace, model_end, "survey_azi_rad", inputs->survey_azi_rad, SOLVER_MAX_SURVEY, &nz) == 0 &&
ni == ns && nz == ns) {
inputs->survey_station_count = ns;
}
}
const char *sc = strstr(buffer, "\"surfaceCard\"");
if (sc && inputs->workflow == 1) {
const char *b = strchr(sc, '{');
if (b) {
const char *end = find_matching_brace(b);
if (end) {
int np = 0, nl = 0;
if (parse_object_array(b, end, "position_m", inputs->surface_position_m, SOLVER_MAX_SURFACE, &np) == 0 &&
parse_object_array(b, end, "load_n", inputs->surface_load_n, SOLVER_MAX_SURFACE, &nl) == 0 && np == nl && np > 0) {
inputs->surface_count = np;
int nt = 0;
if (parse_object_array(b, end, "time_s", inputs->surface_time_s, SOLVER_MAX_SURFACE, &nt) == 0 && nt == np) {
inputs->surface_has_time = 1;
}
}
}
}
}
const char *opt = strstr(buffer, "\"options\"");
if (opt) {
const char *b = strchr(opt, '{');
if (b) {
const char *end = find_matching_brace(b);
if (end) {
(void)parse_object_bool(b, end, "enableProfiles", &inputs->enable_profiles);
(void)parse_object_bool(b, end, "enableDiagnosticsDetail", &inputs->enable_diagnostics_detail);
(void)parse_object_bool(b, end, "enableFourierBaseline", &inputs->enable_fourier_baseline);
(void)parse_object_int(b, end, "fourierHarmonics", &inputs->fourier_harmonics);
}
}
}
if (inputs->fourier_harmonics < 1) inputs->fourier_harmonics = 1;
if (inputs->fourier_harmonics > SOLVER_MAX_FOURIER_HARMONICS) inputs->fourier_harmonics = SOLVER_MAX_FOURIER_HARMONICS;
return 0;
}

177
solver-c/src/main.c Normal file
View File

@@ -0,0 +1,177 @@
#include "solver.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
static void print_json_output(const SolverOutputs *outputs) {
printf("{\n");
printf(" \"pointCount\": %d,\n", outputs->point_count);
printf(" \"maxPolishedLoad\": %.6f,\n", outputs->max_polished_load);
printf(" \"minPolishedLoad\": %.6f,\n", outputs->min_polished_load);
printf(" \"maxDownholeLoad\": %.6f,\n", outputs->max_downhole_load);
printf(" \"minDownholeLoad\": %.6f,\n", outputs->min_downhole_load);
printf(" \"gasInterference\": %s,\n", outputs->gas_interference ? "true" : "false");
printf(" \"maxCfl\": %.6f,\n", outputs->max_cfl);
printf(" \"waveSpeedRefMPerS\": %.6f,\n", outputs->wave_speed_ref_m_s);
printf(" \"warnings\": [");
for (int i = 0; i < outputs->warning_count; i++) {
if (i > 0) printf(", ");
printf("\"%s\"", outputs->warnings[i]);
}
printf("],\n");
printf(" \"card\": [\n");
for (int i = 0; i < outputs->point_count; i++) {
printf(" {\"position\": %.6f, \"polishedLoad\": %.6f, \"downholeLoad\": %.6f, \"polishedStressPa\": %.6f, \"sideLoadN\": %.6f}%s\n",
outputs->position[i],
outputs->polished_load[i],
outputs->downhole_load[i],
outputs->polished_stress_pa[i],
outputs->side_load_profile_n[i],
(i == outputs->point_count - 1) ? "" : ",");
}
printf(" ],\n");
double pmin = outputs->pump_position_m[0];
double pmax = outputs->pump_position_m[0];
for (int i = 1; i < outputs->point_count; i++) {
if (outputs->pump_position_m[i] < pmin) pmin = outputs->pump_position_m[i];
if (outputs->pump_position_m[i] > pmax) pmax = outputs->pump_position_m[i];
}
printf(" \"pumpMovement\": {\n");
printf(" \"stroke\": %.6f,\n", pmax - pmin);
printf(" \"position\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->pump_position_m[i]);
}
printf("],\n");
printf(" \"velocity\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->pump_velocity_m_s[i]);
}
printf("]\n");
printf(" },\n");
printf(" \"profiles\": {\n");
printf(" \"nodeCount\": %d,\n", outputs->profile_node_count);
printf(" \"trajectory3D\": [");
for (int i = 0; i < outputs->profile_node_count; i++) {
if (i > 0) printf(", ");
printf("{\"md\": %.6f, \"curvature\": %.9f, \"inclination\": %.9f, \"azimuth\": %.9f}",
outputs->profile_md_m[i],
outputs->profile_curvature_1pm[i],
outputs->profile_inclination_rad[i],
outputs->profile_azimuth_rad[i]);
}
printf("],\n");
printf(" \"sideLoadProfile\": [");
for (int i = 0; i < outputs->profile_node_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->profile_side_load_n[i]);
}
printf("],\n");
printf(" \"frictionProfile\": [");
for (int i = 0; i < outputs->profile_node_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->profile_friction_n[i]);
}
printf("]\n");
printf(" },\n");
printf(" \"diagnostics\": {\n");
printf(" \"valveStates\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("{\"travelingOpen\": %s, \"standingOpen\": %s}",
outputs->valve_traveling_open[i] ? "true" : "false",
outputs->valve_standing_open[i] ? "true" : "false");
}
printf("],\n");
printf(" \"chamberPressurePa\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->chamber_pressure_pa[i]);
}
printf("],\n");
printf(" \"gasFraction\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->gas_fraction[i]);
}
printf("]\n");
printf(" },\n");
printf(" \"fourierBaseline\": ");
if (outputs->fourier_harmonics_used > 0) {
printf("{\"harmonics\": %d, \"residualRmsPolished\": %.6f, \"residualRmsDownhole\": %.6f, \"card\": [",
outputs->fourier_harmonics_used,
outputs->fourier_residual_rms_polished,
outputs->fourier_residual_rms_downhole);
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("{\"position\": %.6f, \"polishedLoad\": %.6f, \"downholeLoad\": %.6f}",
outputs->position[i],
outputs->fourier_polished_load[i],
outputs->fourier_downhole_load[i]);
}
printf("]}");
} else {
printf("null");
}
printf("\n");
printf("}\n");
}
int main(int argc, char **argv) {
char *buf = NULL;
size_t len = 0;
if (argc > 1 && strcmp(argv[1], "--stdin") == 0) {
if (solver_read_json_stdin(&buf, &len) != 0) {
fprintf(stderr, "failed to read stdin json\n");
return 1;
}
} else if (argc == 2) {
/* allow file path for tests */
FILE *fp = fopen(argv[1], "rb");
if (!fp) {
fprintf(stderr, "usage: solver_main --stdin OR solver_main <path.json>\n");
return 1;
}
fseek(fp, 0, SEEK_END);
long sz = ftell(fp);
fseek(fp, 0, SEEK_SET);
buf = (char *)malloc((size_t)sz + 1);
fread(buf, 1, (size_t)sz, fp);
fclose(fp);
buf[sz] = '\0';
} else {
fprintf(stderr, "usage: solver_main --stdin OR solver_main <path.json>\n");
return 1;
}
SolverInputs inputs;
if (solver_parse_json_inputs(buf, &inputs) != 0) {
fprintf(stderr, "failed to parse json inputs\n");
free(buf);
return 2;
}
free(buf);
SolverOutputs outputs;
int rc = 0;
if (inputs.workflow == 1) {
rc = solver_run_diagnostic_fdm(&inputs, &outputs);
} else {
rc = solver_run_fdm(&inputs, &outputs);
}
if (rc != 0) {
fprintf(stderr, "solver failed: %d\n", rc);
return 3;
}
print_json_output(&outputs);
return 0;
}

170
solver-c/src/main_fea.c Normal file
View File

@@ -0,0 +1,170 @@
#include "solver.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
static void print_json_output(const SolverOutputs *outputs) {
printf("{\n");
printf(" \"pointCount\": %d,\n", outputs->point_count);
printf(" \"maxPolishedLoad\": %.6f,\n", outputs->max_polished_load);
printf(" \"minPolishedLoad\": %.6f,\n", outputs->min_polished_load);
printf(" \"maxDownholeLoad\": %.6f,\n", outputs->max_downhole_load);
printf(" \"minDownholeLoad\": %.6f,\n", outputs->min_downhole_load);
printf(" \"gasInterference\": %s,\n", outputs->gas_interference ? "true" : "false");
printf(" \"maxCfl\": %.6f,\n", outputs->max_cfl);
printf(" \"waveSpeedRefMPerS\": %.6f,\n", outputs->wave_speed_ref_m_s);
printf(" \"warnings\": [");
for (int i = 0; i < outputs->warning_count; i++) {
if (i > 0) printf(", ");
printf("\"%s\"", outputs->warnings[i]);
}
printf("],\n");
printf(" \"card\": [\n");
for (int i = 0; i < outputs->point_count; i++) {
printf(" {\"position\": %.6f, \"polishedLoad\": %.6f, \"downholeLoad\": %.6f, \"polishedStressPa\": %.6f, \"sideLoadN\": %.6f}%s\n",
outputs->position[i],
outputs->polished_load[i],
outputs->downhole_load[i],
outputs->polished_stress_pa[i],
outputs->side_load_profile_n[i],
(i == outputs->point_count - 1) ? "" : ",");
}
printf(" ],\n");
double pmin = outputs->pump_position_m[0];
double pmax = outputs->pump_position_m[0];
for (int i = 1; i < outputs->point_count; i++) {
if (outputs->pump_position_m[i] < pmin) pmin = outputs->pump_position_m[i];
if (outputs->pump_position_m[i] > pmax) pmax = outputs->pump_position_m[i];
}
printf(" \"pumpMovement\": {\n");
printf(" \"stroke\": %.6f,\n", pmax - pmin);
printf(" \"position\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->pump_position_m[i]);
}
printf("],\n");
printf(" \"velocity\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->pump_velocity_m_s[i]);
}
printf("]\n");
printf(" },\n");
printf(" \"profiles\": {\n");
printf(" \"nodeCount\": %d,\n", outputs->profile_node_count);
printf(" \"trajectory3D\": [");
for (int i = 0; i < outputs->profile_node_count; i++) {
if (i > 0) printf(", ");
printf("{\"md\": %.6f, \"curvature\": %.9f, \"inclination\": %.9f, \"azimuth\": %.9f}",
outputs->profile_md_m[i],
outputs->profile_curvature_1pm[i],
outputs->profile_inclination_rad[i],
outputs->profile_azimuth_rad[i]);
}
printf("],\n");
printf(" \"sideLoadProfile\": [");
for (int i = 0; i < outputs->profile_node_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->profile_side_load_n[i]);
}
printf("],\n");
printf(" \"frictionProfile\": [");
for (int i = 0; i < outputs->profile_node_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->profile_friction_n[i]);
}
printf("]\n");
printf(" },\n");
printf(" \"diagnostics\": {\n");
printf(" \"valveStates\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("{\"travelingOpen\": %s, \"standingOpen\": %s}",
outputs->valve_traveling_open[i] ? "true" : "false",
outputs->valve_standing_open[i] ? "true" : "false");
}
printf("],\n");
printf(" \"chamberPressurePa\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->chamber_pressure_pa[i]);
}
printf("],\n");
printf(" \"gasFraction\": [");
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("%.6f", outputs->gas_fraction[i]);
}
printf("]\n");
printf(" },\n");
printf(" \"fourierBaseline\": ");
if (outputs->fourier_harmonics_used > 0) {
printf("{\"harmonics\": %d, \"residualRmsPolished\": %.6f, \"residualRmsDownhole\": %.6f, \"card\": [",
outputs->fourier_harmonics_used,
outputs->fourier_residual_rms_polished,
outputs->fourier_residual_rms_downhole);
for (int i = 0; i < outputs->point_count; i++) {
if (i > 0) printf(", ");
printf("{\"position\": %.6f, \"polishedLoad\": %.6f, \"downholeLoad\": %.6f}",
outputs->position[i],
outputs->fourier_polished_load[i],
outputs->fourier_downhole_load[i]);
}
printf("]}");
} else {
printf("null");
}
printf("\n");
printf("}\n");
}
int main(int argc, char **argv) {
char *buf = NULL;
if (argc > 1 && strcmp(argv[1], "--stdin") == 0) {
size_t len = 0;
if (solver_read_json_stdin(&buf, &len) != 0) {
fprintf(stderr, "failed to read stdin json\n");
return 1;
}
} else if (argc == 2) {
FILE *fp = fopen(argv[1], "rb");
if (!fp) {
fprintf(stderr, "usage: solver_fea_main --stdin OR solver_fea_main <path.json>\n");
return 1;
}
fseek(fp, 0, SEEK_END);
long sz = ftell(fp);
fseek(fp, 0, SEEK_SET);
buf = (char *)malloc((size_t)sz + 1);
fread(buf, 1, (size_t)sz, fp);
fclose(fp);
buf[sz] = '\0';
} else {
fprintf(stderr, "usage: solver_fea_main --stdin OR solver_fea_main <path.json>\n");
return 1;
}
SolverInputs inputs;
if (solver_parse_json_inputs(buf, &inputs) != 0) {
fprintf(stderr, "failed to parse json inputs\n");
free(buf);
return 2;
}
free(buf);
SolverOutputs outputs;
if (solver_run_fea(&inputs, &outputs) != 0) {
fprintf(stderr, "solver failed\n");
return 3;
}
print_json_output(&outputs);
return 0;
}

207
solver-c/src/solver.c Normal file
View File

@@ -0,0 +1,207 @@
#include "solver.h"
#include "solver_internal.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define FDM_NX 48
static void build_default_rod_nodes(SolverInputs *inputs, int nx, double rod_length_m) {
const int nodes = nx + 1;
const double E = 2.05e11;
const double rho = 7850.0;
const double taper = solver_input_or_default(inputs->taper_factor, 1.0);
const double area = solver_clamp((0.00042 + inputs->rod_friction_coefficient * 0.00008) * taper, 0.00025, 0.0009);
inputs->has_variable_rod = 1;
inputs->rod_node_count = nodes;
for (int i = 0; i < nodes; i++) {
inputs->area_m2[i] = area;
inputs->modulus_pa[i] = E;
inputs->density_kg_m3[i] = rho;
}
(void)rod_length_m;
}
static double segment_mean(const double *a, int i, int max_i) {
const int j = i + 1 > max_i ? max_i : i + 1;
return 0.5 * (a[i] + a[j]);
}
/*
* Explicit FDM core matches legacy stable implementation (pre-JSON refactor).
* Variable rod arrays are populated for API/verbose; dynamics use scalar EA.
*/
int solver_run_fdm(const SolverInputs *inputs_in, SolverOutputs *outputs) {
if (inputs_in == NULL || outputs == NULL) {
return -1;
}
SolverInputs inputs_local = *inputs_in;
SolverInputs *inputs = &inputs_local;
memset(outputs, 0, sizeof(SolverOutputs));
outputs->point_count = 200;
solver_add_warning(outputs, "FDM model: damped 1D wave equation solved in time-space grid");
if (inputs->pumping_speed <= 0.0) {
solver_add_warning(outputs, "Non-positive pumping speed, using fallback 1.0 SPM");
}
if (inputs->pump_depth <= 0.0) {
solver_add_warning(outputs, "Non-positive pump depth, using fallback depth");
}
const double spm = solver_input_or_default(inputs->pumping_speed, 1.0);
const double depth = solver_input_or_default(inputs->pump_depth, 1000.0);
const double anchor = solver_input_or_default(inputs->tubing_anchor_location, depth * 0.8);
const double rod_length = solver_clamp(depth - anchor, 250.0, 3500.0);
const double period = 60.0 / spm;
const double dt = period / (double)(outputs->point_count - 1);
const double dx = rod_length / (double)FDM_NX;
const int nx = FDM_NX;
if (!inputs->has_variable_rod || inputs->rod_node_count != nx + 1) {
build_default_rod_nodes(inputs, nx, rod_length);
}
solver_trajectory_preprocess(inputs, nx, rod_length);
const double E = 2.05e11;
const double rho = 7850.0;
const double taper_factor = solver_clamp(solver_input_or_default(inputs->taper_factor, 1.0), 0.65, 1.25);
const double traj_fric_mul = solver_clamp(solver_input_or_default(inputs->trajectory_friction_multiplier, 1.0), 1.0, 1.8);
const double area = solver_clamp((0.00042 + inputs->rod_friction_coefficient * 0.00008) * taper_factor, 0.00025, 0.0009);
const double a = sqrt(E / rho);
outputs->wave_speed_ref_m_s = a;
const double cfl = a * dt / dx;
/* Tighter explicit stability margin than legacy 0.95 for stiff pump BC + coarse grid */
const double cfl_eff = (cfl > 0.95) ? 0.95 : cfl;
outputs->max_cfl = cfl;
const double r = cfl_eff * cfl_eff;
const double gamma = solver_clamp(0.08 + inputs->rod_friction_coefficient * 0.25, 0.05, 0.3);
const double damping = gamma * dt;
const double surf_amp = solver_clamp(depth * 0.00045, 0.4, 2.6);
const double baseline = depth * 2.05 + 3900.0;
const double k_pump = solver_clamp(1.8e5 + depth * 180.0, 1.2e5, 9.0e5);
const double c_pump = solver_clamp(220.0 + inputs->rod_friction_coefficient * 800.0, 120.0, 1200.0);
const double pump_friction = solver_clamp(inputs->pump_friction * traj_fric_mul, 0.0, 7000.0);
const double stuffing_friction = solver_clamp(inputs->stuffing_box_friction * traj_fric_mul, 0.0, 3500.0);
double side_nodes[SOLVER_MAX_NODES];
double fric_nodes[SOLVER_MAX_NODES];
memset(side_nodes, 0, sizeof(side_nodes));
memset(fric_nodes, 0, sizeof(fric_nodes));
if (spm > 20.0 || cfl > 1.0) {
solver_add_warning(outputs, "FDM settings near CFL limit; solution uses stabilized time step");
}
if (fabs(taper_factor - 1.0) > 1e-6 || fabs(traj_fric_mul - 1.0) > 1e-6) {
solver_add_warning(outputs, "FDM using taper and trajectory-coupled coefficients");
}
solver_init_output_ranges(outputs);
double u_prev[FDM_NX + 1];
double u_curr[FDM_NX + 1];
double u_next[FDM_NX + 1];
memset(u_prev, 0, sizeof(u_prev));
memset(u_curr, 0, sizeof(u_curr));
memset(u_next, 0, sizeof(u_next));
const double omega = 2.0 * M_PI / period;
for (int n = 0; n < outputs->point_count; n++) {
const double t = n * dt;
const double u0 = surf_amp * sin(omega * t);
const double v0 = surf_amp * omega * cos(omega * t);
u_curr[0] = u0;
if (n == 0) {
for (int i = 1; i <= nx; i++) {
u_curr[i] = u0 * (1.0 - (double)i / (double)nx);
u_prev[i] = u_curr[i];
}
}
for (int i = 1; i < nx; i++) {
const double lap = u_curr[i + 1] - 2.0 * u_curr[i] + u_curr[i - 1];
const double v_i = (u_curr[i] - u_prev[i]) / dt;
const double t_i = E * area * (u_curr[i] - u_curr[i - 1]) / dx;
const double side_i = solver_compute_side_load_node(inputs, t_i, i, dx);
const double fric_i = solver_compute_friction_node(inputs, side_i, v_i, i);
side_nodes[i] = side_i;
fric_nodes[i] = fric_i;
const double body = -(side_i + fric_i) * dx / fmax(E * area, 1e-6);
u_next[i] = (2.0 - damping) * u_curr[i] - (1.0 - damping) * u_prev[i] + r * lap + body;
}
const double uN = u_curr[nx];
const double vN = (u_curr[nx] - u_prev[nx]) / dt;
const double tension_end = E * area * (u_curr[nx] - u_curr[nx - 1]) / dx;
const double side_mid = solver_compute_side_load_node(inputs, tension_end, nx, dx);
const double friction_mid = solver_compute_friction_node(inputs, side_mid, vN, nx);
side_nodes[nx] = side_mid;
fric_nodes[nx] = friction_mid;
const double spring = k_pump * uN;
const double visc = c_pump * vN;
const double fric = (pump_friction + 0.65 * stuffing_friction) * solver_signum(vN) + friction_mid;
const double pump_force = spring + visc + fric - side_mid;
u_next[nx] = u_next[nx - 1] + (pump_force * dx) / (E * area);
u_next[0] = surf_amp * sin(omega * (t + dt));
const double e01 = segment_mean(inputs->modulus_pa, 0, nx);
const double a01 = segment_mean(inputs->area_m2, 0, nx);
const double polished_tension = E * area * (u_curr[1] - u_curr[0]) / dx;
const double polished = baseline + polished_tension + stuffing_friction * solver_signum(v0);
const double downhole = baseline + pump_force;
outputs->position[n] = u_curr[0];
outputs->polished_load[n] = polished;
outputs->downhole_load[n] = downhole;
outputs->pump_position_m[n] = u_curr[nx];
outputs->polished_stress_pa[n] = polished_tension / fmax(a01, 1e-12);
const double side_top = solver_compute_side_load_node(inputs, polished_tension, 0, dx);
side_nodes[0] = side_top;
fric_nodes[0] = solver_compute_friction_node(inputs, side_top, v0, 0);
outputs->side_load_profile_n[n] = side_top;
(void)e01;
solver_update_output_ranges(outputs, polished, downhole);
for (int i = 0; i <= nx; i++) {
u_prev[i] = u_curr[i];
u_curr[i] = u_next[i];
}
}
outputs->gas_interference =
(inputs->pump_fillage_option == 2 || (inputs->percent_pump_fillage > 0.0 && inputs->percent_pump_fillage < 88.0))
? 1
: 0;
for (int i = 0; i < outputs->point_count; i++) {
if (i == 0) {
outputs->pump_velocity_m_s[i] = 0.0;
} else {
outputs->pump_velocity_m_s[i] = (outputs->pump_position_m[i] - outputs->pump_position_m[i - 1]) / dt;
}
}
if (outputs->point_count > 1) {
outputs->pump_velocity_m_s[0] = outputs->pump_velocity_m_s[1];
}
for (int i = 0; i < outputs->point_count; i++) {
solver_valve_state_step(inputs, outputs, i, outputs->pump_position_m[i], outputs->pump_velocity_m_s[i], outputs->downhole_load[i]);
}
solver_fill_profiles(inputs, outputs, nx + 1, rod_length, side_nodes, fric_nodes);
if (inputs->enable_fourier_baseline) {
(void)solver_compute_fourier_baseline(inputs, outputs);
}
return 0;
}
int solver_run(const SolverInputs *inputs, SolverOutputs *outputs) {
return solver_run_fdm(inputs, outputs);
}

117
solver-c/src/solver_common.c Normal file
View File

@@ -0,0 +1,117 @@
#include "solver_internal.h"
#include <math.h>
#include <string.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
double solver_clamp(double v, double lo, double hi) {
if (v < lo) return lo;
if (v > hi) return hi;
return v;
}
double solver_signum(double v) {
if (v > 0.0) return 1.0;
if (v < 0.0) return -1.0;
return 0.0;
}
void solver_add_warning(SolverOutputs *outputs, const char *msg) {
if (outputs->warning_count >= SOLVER_MAX_WARNINGS) {
return;
}
strncpy(outputs->warnings[outputs->warning_count], msg, SOLVER_WARNING_LEN - 1);
outputs->warnings[outputs->warning_count][SOLVER_WARNING_LEN - 1] = '\0';
outputs->warning_count += 1;
}
void solver_init_output_ranges(SolverOutputs *outputs) {
outputs->max_polished_load = -1e99;
outputs->min_polished_load = 1e99;
outputs->max_downhole_load = -1e99;
outputs->min_downhole_load = 1e99;
}
void solver_update_output_ranges(SolverOutputs *outputs, double polished, double downhole) {
if (polished > outputs->max_polished_load) outputs->max_polished_load = polished;
if (polished < outputs->min_polished_load) outputs->min_polished_load = polished;
if (downhole > outputs->max_downhole_load) outputs->max_downhole_load = downhole;
if (downhole < outputs->min_downhole_load) outputs->min_downhole_load = downhole;
}
double solver_input_or_default(double value, double fallback) {
return (value > 0.0) ? value : fallback;
}
double solver_compute_side_load_node(const SolverInputs *inputs, double tension_n, int node_idx, double ds) {
if (!inputs) return 0.0;
const int i = (node_idx < 0) ? 0 : ((node_idx >= SOLVER_MAX_NODES) ? SOLVER_MAX_NODES - 1 : node_idx);
const double kappa = fmax(inputs->node_curvature[i], 0.0);
const double inc = fabs(inputs->node_inc_rad[i]);
const double rho = fmax(inputs->density_kg_m3[i], 1.0);
const double area = fmax(inputs->area_m2[i], 1e-10);
double buoyed_weight = fmax((rho - inputs->fluid_density_kg_m3) * area * inputs->gravity, 0.0);
if (inputs->weight_n_per_m[i] > 0.0) {
buoyed_weight = fmax(inputs->weight_n_per_m[i] - inputs->fluid_density_kg_m3 * area * inputs->gravity, 0.0);
}
const double guide_weight = fmax(inputs->guide_weight_n_per_m[i], 0.0);
const double sinker_area = M_PI * 0.25 * inputs->sinker_bar_diameter_m * inputs->sinker_bar_diameter_m;
const double sinker_weight = (inputs->sinker_bar_length_m > 0.0)
? fmax((7850.0 - inputs->fluid_density_kg_m3) * sinker_area * inputs->gravity, 0.0)
: 0.0;
/* Lukasiewicz-inspired normal force combination: curvature tension + gravity/inclination lateral component */
const double normal = fabs(tension_n) * kappa + (buoyed_weight + guide_weight + sinker_weight) * sin(inc) * fmax(ds, 1e-6);
return fmax(normal, 0.0);
}
double solver_compute_friction_node(const SolverInputs *inputs, double side_load_n, double velocity_m_s, int node_idx) {
if (!inputs) return 0.0;
const int i = (node_idx < 0) ? 0 : ((node_idx >= SOLVER_MAX_NODES) ? SOLVER_MAX_NODES - 1 : node_idx);
const double base_mu = solver_clamp(inputs->rod_friction_coefficient, 0.0, 1.5);
const double guide_scale = (inputs->molded_guide_mu_scale + inputs->wheeled_guide_mu_scale + inputs->other_guide_mu_scale) / 3.0;
const double inc_scale = 1.0 + 0.5 * fabs(sin(inputs->node_inc_rad[i]));
const double mts_scale = inputs->mts_n[i] > 0.0 ? solver_clamp(inputs->mts_n[i] / 8.0e5, 0.5, 1.5) : 1.0;
const double mu = base_mu * solver_clamp(guide_scale, 0.1, 3.0) * inc_scale * mts_scale;
return mu * fmax(side_load_n, 0.0) * solver_signum(velocity_m_s);
}
void solver_fill_profiles(const SolverInputs *inputs, SolverOutputs *outputs, int node_count, double rod_length_m,
const double *side_load_nodes, const double *friction_nodes) {
if (!inputs || !outputs) return;
const int n = (node_count > SOLVER_MAX_NODES) ? SOLVER_MAX_NODES : node_count;
outputs->profile_node_count = n;
for (int i = 0; i < n; i++) {
outputs->profile_md_m[i] = (n > 1) ? rod_length_m * (double)i / (double)(n - 1) : 0.0;
outputs->profile_curvature_1pm[i] = inputs->node_curvature[i];
outputs->profile_inclination_rad[i] = inputs->node_inc_rad[i];
outputs->profile_azimuth_rad[i] = inputs->node_azi_rad[i];
outputs->profile_side_load_n[i] = side_load_nodes ? side_load_nodes[i] : 0.0;
outputs->profile_friction_n[i] = friction_nodes ? friction_nodes[i] : 0.0;
}
}
void solver_valve_state_step(const SolverInputs *inputs, SolverOutputs *outputs, int step_idx, double pump_position_m,
double pump_velocity_m_s, double downhole_load_n) {
if (!inputs || !outputs) return;
if (step_idx < 0 || step_idx >= SOLVER_MAX_POINTS) return;
const double area = M_PI * fmax(inputs->pump_diameter_m, 0.02) * fmax(inputs->pump_diameter_m, 0.02) * 0.25;
const double hydro = inputs->pump_intake_pressure_pa + inputs->fluid_density_kg_m3 * inputs->gravity * fmax(inputs->pump_depth, 0.0);
const double piston_pressure = hydro + downhole_load_n / fmax(area, 1e-6);
const double discharge_pressure = hydro + 0.5 * fabs(downhole_load_n) / fmax(area, 1e-6);
const int standing_open = (piston_pressure < hydro * 1.02) ? 1 : 0;
const int traveling_open = (pump_velocity_m_s < 0.0 && piston_pressure > discharge_pressure * 0.98) ? 1 : 0;
outputs->valve_standing_open[step_idx] = standing_open;
outputs->valve_traveling_open[step_idx] = traveling_open;
const double stroke_norm = solver_clamp(fabs(pump_position_m) / fmax(fabs(pump_position_m) + 1e-6, 1e-6), 0.0, 1.0);
const double gas_base = solver_clamp(1.0 - inputs->percent_pump_fillage / 100.0, 0.0, 0.9);
double gas = gas_base + 0.2 * (1.0 - stroke_norm) + ((standing_open && traveling_open) ? 0.05 : 0.0);
gas = solver_clamp(gas, 0.0, 0.98);
outputs->gas_fraction[step_idx] = gas;
outputs->chamber_pressure_pa[step_idx] = piston_pressure * (1.0 + 0.25 * gas);
if (gas > 0.35) outputs->gas_interference = 1;
}

239
solver-c/src/solver_diagnostic.c Normal file
View File

@@ -0,0 +1,239 @@
#include "solver.h"
#include "solver_internal.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define FDM_NX 48
static double mean_series(const double *v, int n) {
if (n <= 0) {
return 0.0;
}
double s = 0.0;
for (int i = 0; i < n; i++) {
s += v[i];
}
return s / (double)n;
}
static double segment_mean(const double *a, int i, int max_i) {
const int j = i + 1 > max_i ? max_i : i + 1;
return 0.5 * (a[i] + a[j]);
}
int solver_run_diagnostic_fdm(const SolverInputs *inputs_in, SolverOutputs *outputs) {
if (inputs_in == NULL || outputs == NULL) {
return -1;
}
SolverInputs inputs_local = *inputs_in;
SolverInputs *inputs = &inputs_local;
memset(outputs, 0, sizeof(SolverOutputs));
solver_add_warning(outputs, "Diagnostic FDM: surface card BC (Everitt & Jennings stencil, MATH_SPEC.md)");
const int point_count = inputs->surface_count;
if (point_count < 40) {
return -2;
}
outputs->point_count = point_count;
const double pumping_speed = inputs->pumping_speed > 0.0 ? inputs->pumping_speed : 5.0;
const double period = 60.0 / pumping_speed;
double dt = period / (double)(point_count - 1);
if (inputs->surface_has_time && point_count > 1) {
dt = (inputs->surface_time_s[point_count - 1] - inputs->surface_time_s[0]) / (double)(point_count - 1);
if (dt <= 0.0 || dt > period * 2.0) {
dt = period / (double)(point_count - 1);
}
}
const double depth = inputs->pump_depth > 0.0 ? inputs->pump_depth : 1000.0;
const double anchor = inputs->tubing_anchor_location > 0.0 ? inputs->tubing_anchor_location : depth * 0.8;
const double rod_length = solver_clamp(depth - anchor, 250.0, 3500.0);
const int nx = FDM_NX;
const int nodes = nx + 1;
const double dx = rod_length / (double)nx;
if (!inputs->has_variable_rod || inputs->rod_node_count != nodes) {
inputs->has_variable_rod = 1;
inputs->rod_node_count = nodes;
const double E = 2.05e11;
const double rho = 7850.0;
const double taper = solver_input_or_default(inputs->taper_factor, 1.0);
const double area = solver_clamp((0.00042 + inputs->rod_friction_coefficient * 0.00008) * taper, 0.00025, 0.0009);
for (int i = 0; i < nodes; i++) {
inputs->area_m2[i] = area;
inputs->modulus_pa[i] = E;
inputs->density_kg_m3[i] = rho;
}
}
solver_trajectory_preprocess(inputs, nx, rod_length);
double e_ref = 0.0;
double rho_ref = 0.0;
for (int i = 0; i < nodes; i++) {
e_ref += inputs->modulus_pa[i];
rho_ref += inputs->density_kg_m3[i];
}
e_ref /= (double)nodes;
rho_ref /= (double)nodes;
const double a = sqrt(e_ref / fmax(rho_ref, 1e-6));
const double cfl = (a * dt) / dx;
const double cfl_eff = (cfl > 0.95) ? 0.95 : cfl;
outputs->max_cfl = cfl;
outputs->wave_speed_ref_m_s = a;
const double r = cfl_eff * cfl_eff;
const double gamma = solver_clamp(0.08 + inputs->rod_friction_coefficient * 0.25, 0.05, 0.3);
const double damping = gamma * dt;
double baseline = depth * 2.05 + 3900.0;
for (int i = 0; i < nx; i++) {
const double rho_i = inputs->density_kg_m3[i];
const double A_i = segment_mean(inputs->area_m2, i, nodes - 1);
const double inc = 0.5 * (inputs->node_inc_rad[i] + inputs->node_inc_rad[i + 1]);
baseline += 0.5 * (rho_i - inputs->fluid_density_kg_m3) * A_i * inputs->gravity * cos(inc) * dx;
}
const double k_pump = solver_clamp(1.8e5 + depth * 180.0, 1.2e5, 9.0e5);
const double c_pump = solver_clamp(220.0 + inputs->rod_friction_coefficient * 800.0, 120.0, 1200.0);
double kappa_mean = 0.0;
for (int i = 0; i < nodes; i++) {
kappa_mean += inputs->node_curvature[i];
}
kappa_mean /= (double)nodes;
const double mu_geo =
inputs->rod_friction_coefficient *
(inputs->molded_guide_mu_scale + inputs->wheeled_guide_mu_scale + inputs->other_guide_mu_scale) / 3.0;
double friction_geom = 1.0 + solver_clamp(mu_geo * kappa_mean * rod_length * 0.08, 0.0, 1.2);
if (!inputs->geometry_valid) {
friction_geom = solver_clamp(inputs->trajectory_friction_multiplier, 1.0, 1.8);
}
const double pump_friction = solver_clamp(inputs->pump_friction * friction_geom, 0.0, 6000.0);
const double stuffing_friction = solver_clamp(inputs->stuffing_box_friction * friction_geom, 0.0, 2800.0);
double side_nodes[SOLVER_MAX_NODES];
double fric_nodes[SOLVER_MAX_NODES];
memset(side_nodes, 0, sizeof(side_nodes));
memset(fric_nodes, 0, sizeof(fric_nodes));
double u_prev[SOLVER_MAX_NODES];
double u_curr[SOLVER_MAX_NODES];
double u_next[SOLVER_MAX_NODES];
memset(u_prev, 0, sizeof(u_prev));
memset(u_curr, 0, sizeof(u_curr));
memset(u_next, 0, sizeof(u_next));
const double pos_mean = mean_series(inputs->surface_position_m, point_count);
const double load_mean = mean_series(inputs->surface_load_n, point_count);
solver_init_output_ranges(outputs);
for (int n = 0; n < point_count; n++) {
const double pos = inputs->surface_position_m[n] - pos_mean;
const double load = inputs->surface_load_n[n] - load_mean;
u_curr[0] = pos;
const double e01 = segment_mean(inputs->modulus_pa, 0, nodes - 1);
const double a01 = segment_mean(inputs->area_m2, 0, nodes - 1);
u_curr[1] = u_curr[0] + (load * dx) / (e01 * a01 + 1e-30);
if (n == 0) {
for (int i = 2; i <= nx; i++) {
const double ratio = 1.0 - (double)i / (double)nx;
u_curr[i] = u_curr[1] * ratio;
u_prev[i] = u_curr[i];
}
u_prev[0] = u_curr[0];
u_prev[1] = u_curr[1];
}
for (int i = 1; i < nx; i++) {
const double lap = u_curr[i + 1] - 2.0 * u_curr[i] + u_curr[i - 1];
const double v_i = (u_curr[i] - u_prev[i]) / dt;
const double e_i = segment_mean(inputs->modulus_pa, i - 1, nodes - 1);
const double a_i = segment_mean(inputs->area_m2, i - 1, nodes - 1);
const double t_i = e_i * a_i * (u_curr[i] - u_curr[i - 1]) / dx;
const double side_i = solver_compute_side_load_node(inputs, t_i, i, dx);
const double fric_i = solver_compute_friction_node(inputs, side_i, v_i, i);
side_nodes[i] = side_i;
fric_nodes[i] = fric_i;
const double body = -(side_i + fric_i) * dx / fmax(e_i * a_i, 1e-6);
u_next[i] = (2.0 - damping) * u_curr[i] - (1.0 - damping) * u_prev[i] + r * lap + body;
}
const double e_end = segment_mean(inputs->modulus_pa, nx - 1, nodes - 1);
const double a_end = segment_mean(inputs->area_m2, nx - 1, nodes - 1);
const double vN = (u_curr[nx] - u_prev[nx]) / dt;
const double tension_end = e_end * a_end * (u_curr[nx] - u_curr[nx - 1]) / dx;
const double side_mid = solver_compute_side_load_node(inputs, tension_end, nx, dx);
const double fric_mid = solver_compute_friction_node(inputs, side_mid, vN, nx);
side_nodes[nx] = side_mid;
fric_nodes[nx] = fric_mid;
const double spring = k_pump * u_curr[nx];
const double visc = c_pump * vN;
const double fric = (pump_friction + 0.65 * stuffing_friction) * solver_signum(vN) + fric_mid;
const double pump_force = spring + visc + fric - side_mid;
u_next[nx] = u_next[nx - 1] + (pump_force * dx) / (e_end * a_end + 1e-30);
const int next_idx = (n + 1) % point_count;
const double pos_next = inputs->surface_position_m[next_idx] - pos_mean;
const double load_next = inputs->surface_load_n[next_idx] - load_mean;
u_next[0] = pos_next;
u_next[1] = u_next[0] + (load_next * dx) / (e01 * a01 + 1e-30);
const double v0 = (u_curr[0] - u_prev[0]) / dt;
const double t_surf = (e01 * a01 * (u_curr[1] - u_curr[0])) / dx;
const double polished = baseline + t_surf + stuffing_friction * solver_signum(v0);
const double downhole = baseline + pump_force;
outputs->position[n] = u_curr[0];
outputs->polished_load[n] = polished;
outputs->downhole_load[n] = downhole;
outputs->pump_position_m[n] = u_curr[nx];
outputs->polished_stress_pa[n] = t_surf / fmax(a01, 1e-12);
const double side_top = solver_compute_side_load_node(inputs, t_surf, 0, dx);
side_nodes[0] = side_top;
fric_nodes[0] = solver_compute_friction_node(inputs, side_top, v0, 0);
outputs->side_load_profile_n[n] = side_top;
solver_update_output_ranges(outputs, polished, downhole);
for (int i = 0; i <= nx; i++) {
u_prev[i] = u_curr[i];
u_curr[i] = u_next[i];
}
}
outputs->gas_interference =
(inputs->pump_fillage_option == 2 || (inputs->percent_pump_fillage > 0.0 && inputs->percent_pump_fillage < 88.0))
? 1
: 0;
for (int i = 0; i < point_count; i++) {
if (i == 0) {
outputs->pump_velocity_m_s[i] = 0.0;
} else {
outputs->pump_velocity_m_s[i] = (outputs->pump_position_m[i] - outputs->pump_position_m[i - 1]) / dt;
}
}
if (point_count > 1) {
outputs->pump_velocity_m_s[0] = outputs->pump_velocity_m_s[1];
}
for (int i = 0; i < point_count; i++) {
solver_valve_state_step(inputs, outputs, i, outputs->pump_position_m[i], outputs->pump_velocity_m_s[i], outputs->downhole_load[i]);
}
solver_fill_profiles(inputs, outputs, nodes, rod_length, side_nodes, fric_nodes);
if (inputs->enable_fourier_baseline) {
(void)solver_compute_fourier_baseline(inputs, outputs);
}
return 0;
}

432
solver-c/src/solver_fea.c Normal file
View File

@@ -0,0 +1,432 @@
#include "solver.h"
#include "solver_internal.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define FEA_ELEMENTS 24
#define MAX_FEA_NODES 65
static int solve_linear_system(int n, double A[MAX_FEA_NODES][MAX_FEA_NODES], double b[MAX_FEA_NODES], double x[MAX_FEA_NODES]) {
for (int i = 0; i < n; i++) {
int pivot = i;
double max_val = fabs(A[i][i]);
for (int k = i + 1; k < n; k++) {
double val = fabs(A[k][i]);
if (val > max_val) {
max_val = val;
pivot = k;
}
}
if (max_val < 1e-18) {
return -1;
}
if (pivot != i) {
for (int j = i; j < n; j++) {
double tmp = A[i][j];
A[i][j] = A[pivot][j];
A[pivot][j] = tmp;
}
double bt = b[i];
b[i] = b[pivot];
b[pivot] = bt;
}
double diag = A[i][i];
for (int j = i; j < n; j++) {
A[i][j] /= diag;
}
b[i] /= diag;
for (int k = i + 1; k < n; k++) {
double factor = A[k][i];
for (int j = i; j < n; j++) {
A[k][j] -= factor * A[i][j];
}
b[k] -= factor * b[i];
}
}
for (int i = n - 1; i >= 0; i--) {
double sum = b[i];
for (int j = i + 1; j < n; j++) {
sum -= A[i][j] * x[j];
}
x[i] = sum;
}
return 0;
}
static void fea_assemble_system(const SolverInputs *inputs, int nodes, int elements, double dx, double K[MAX_FEA_NODES][MAX_FEA_NODES],
double M[MAX_FEA_NODES][MAX_FEA_NODES], double *alpha_rayleigh, double *beta_rayleigh) {
memset(K, 0, sizeof(double) * MAX_FEA_NODES * MAX_FEA_NODES);
memset(M, 0, sizeof(double) * MAX_FEA_NODES * MAX_FEA_NODES);
double rho_bar = 0.0;
double e_bar = 0.0;
for (int i = 0; i < nodes; i++) {
rho_bar += inputs->density_kg_m3[i];
e_bar += inputs->modulus_pa[i];
}
rho_bar /= (double)nodes;
e_bar /= (double)nodes;
for (int e = 0; e < elements; e++) {
const int i = e;
const int j = e + 1;
const double Ee = 0.5 * (inputs->modulus_pa[i] + inputs->modulus_pa[j]);
const double Ae = 0.5 * (inputs->area_m2[i] + inputs->area_m2[j]);
const double rhoe = 0.5 * (inputs->density_kg_m3[i] + inputs->density_kg_m3[j]);
const double ke = Ee * Ae / dx;
const double me = rhoe * Ae * dx / 6.0;
K[i][i] += ke;
K[i][j] -= ke;
K[j][i] -= ke;
K[j][j] += ke;
M[i][i] += 2.0 * me;
M[i][j] += me;
M[j][i] += me;
M[j][j] += 2.0 * me;
}
const double L = dx * (double)elements;
const double nu = solver_clamp(inputs->non_dim_damping, 0.05, 1.5);
const double omega1 = M_PI * sqrt(e_bar / fmax(rho_bar, 1e-9)) / fmax(L, 1.0);
const double zeta = solver_clamp(0.02 + inputs->upstroke_damping * 0.05 + nu * 0.03, 0.01, 0.25);
*alpha_rayleigh = 2.0 * zeta * omega1;
*beta_rayleigh = zeta / (2.0 * fmax(omega1, 1e-6));
}
static int fea_newmark_step(const SolverInputs *inputs, int nodes, int elements, double dx, double dt, double k_pump, double c_pump,
double pump_friction, double stuffing_friction, double u0, double v0, double bottom_force_extra, double *u,
double *v, double *a_vec, double *polished_tension_out, double *downhole_tension_out) {
(void)elements;
double K[MAX_FEA_NODES][MAX_FEA_NODES];
double M[MAX_FEA_NODES][MAX_FEA_NODES];
double C[MAX_FEA_NODES][MAX_FEA_NODES];
double alpha_r = 0.0;
double beta_r = 0.0;
fea_assemble_system(inputs, nodes, elements, dx, K, M, &alpha_r, &beta_r);
for (int i = 0; i < nodes; i++) {
for (int j = 0; j < nodes; j++) {
C[i][j] = alpha_r * M[i][j] + beta_r * K[i][j];
}
}
K[nodes - 1][nodes - 1] += k_pump;
C[nodes - 1][nodes - 1] += c_pump;
const double beta_nm = 0.25;
const double gamma_nm = 0.5;
const double a0 = 1.0 / (beta_nm * dt * dt);
const double a1 = gamma_nm / (beta_nm * dt);
const double a2 = 1.0 / (beta_nm * dt);
const double a3 = (1.0 / (2.0 * beta_nm)) - 1.0;
const double a4 = (gamma_nm / beta_nm) - 1.0;
const double a5 = dt * ((gamma_nm / (2.0 * beta_nm)) - 1.0);
const double friction_force = (pump_friction + 0.65 * stuffing_friction) * solver_signum(v[nodes - 1]);
double F[MAX_FEA_NODES];
memset(F, 0, sizeof(F));
for (int i = 1; i < nodes; i++) {
const double tension_i = (i > 0) ? (0.5 * (inputs->modulus_pa[i] + inputs->modulus_pa[i - 1]) *
0.5 * (inputs->area_m2[i] + inputs->area_m2[i - 1]) *
(u[i] - u[i - 1]) / fmax(dx, 1e-9))
: 0.0;
const double side_i = solver_compute_side_load_node(inputs, tension_i, i, dx);
const double fric_i = solver_compute_friction_node(inputs, side_i, v[i], i);
F[i] -= side_i + fric_i;
}
F[nodes - 1] -= friction_force + bottom_force_extra;
const int free_nodes = nodes - 1;
double Kff[MAX_FEA_NODES][MAX_FEA_NODES];
double Mff[MAX_FEA_NODES][MAX_FEA_NODES];
double Cff[MAX_FEA_NODES][MAX_FEA_NODES];
for (int i = 0; i < free_nodes; i++) {
for (int j = 0; j < free_nodes; j++) {
Kff[i][j] = K[i + 1][j + 1];
Mff[i][j] = M[i + 1][j + 1];
Cff[i][j] = C[i + 1][j + 1];
}
}
double rhs[MAX_FEA_NODES];
double Keff[MAX_FEA_NODES][MAX_FEA_NODES];
memset(rhs, 0, sizeof(rhs));
memset(Keff, 0, sizeof(Keff));
for (int i = 0; i < free_nodes; i++) {
const int gi = i + 1;
rhs[i] = F[gi] - K[gi][0] * u0 - C[gi][0] * v0;
for (int j = 0; j < free_nodes; j++) {
Keff[i][j] = Kff[i][j] + a0 * Mff[i][j] + a1 * Cff[i][j];
rhs[i] += Mff[i][j] * (a0 * u[j + 1] + a2 * v[j + 1] + a3 * a_vec[j + 1]);
rhs[i] += Cff[i][j] * (a1 * u[j + 1] + a4 * v[j + 1] + a5 * a_vec[j + 1]);
}
}
double u_new[MAX_FEA_NODES];
memset(u_new, 0, sizeof(u_new));
if (solve_linear_system(free_nodes, Keff, rhs, u_new) != 0) {
return -2;
}
double a_new[MAX_FEA_NODES];
double v_new[MAX_FEA_NODES];
memset(a_new, 0, sizeof(a_new));
memset(v_new, 0, sizeof(v_new));
double old_u[MAX_FEA_NODES];
double old_v[MAX_FEA_NODES];
double old_a[MAX_FEA_NODES];
memcpy(old_u, u, sizeof(old_u));
memcpy(old_v, v, sizeof(old_v));
memcpy(old_a, a_vec, sizeof(old_a));
u[0] = u0;
v[0] = v0;
for (int i = 0; i < free_nodes; i++) {
const int gi = i + 1;
u[gi] = u_new[i];
a_new[gi] = a0 * (u[gi] - old_u[gi]) - a2 * old_v[gi] - a3 * old_a[gi];
v_new[gi] = old_v[gi] + dt * ((1.0 - gamma_nm) * old_a[gi] + gamma_nm * a_new[gi]);
}
for (int gi = 1; gi < nodes; gi++) {
a_vec[gi] = a_new[gi];
v[gi] = v_new[gi];
}
const double e01 = 0.5 * (inputs->modulus_pa[0] + inputs->modulus_pa[1]);
const double a01 = 0.5 * (inputs->area_m2[0] + inputs->area_m2[1]);
const double polished_tension = e01 * a01 * (u[1] - u[0]) / dx;
const double e_end = 0.5 * (inputs->modulus_pa[nodes - 1] + inputs->modulus_pa[nodes - 2]);
const double a_end = 0.5 * (inputs->area_m2[nodes - 1] + inputs->area_m2[nodes - 2]);
const double downhole_tension =
e_end * a_end * (u[nodes - 1] - u[nodes - 2]) / dx + k_pump * u[nodes - 1] + c_pump * v[nodes - 1] + friction_force + bottom_force_extra;
*polished_tension_out = polished_tension;
*downhole_tension_out = downhole_tension;
return 0;
}
int solver_run_fea(const SolverInputs *inputs_in, SolverOutputs *outputs) {
if (inputs_in == NULL || outputs == NULL) {
return -1;
}
SolverInputs inputs_local = *inputs_in;
SolverInputs *inputs = &inputs_local;
memset(outputs, 0, sizeof(SolverOutputs));
outputs->point_count = 200;
solver_add_warning(outputs, "FEA: dynamic 1D bar, Newmark-beta, Rayleigh damping, variable EA (Eisner-style bar FEM)");
const double spm = solver_input_or_default(inputs->pumping_speed, 1.0);
const double depth = solver_input_or_default(inputs->pump_depth, 1000.0);
const double anchor = solver_input_or_default(inputs->tubing_anchor_location, depth * 0.8);
const double rod_length = solver_clamp(depth - anchor, 250.0, 3500.0);
const int elements = FEA_ELEMENTS;
const int nodes = elements + 1;
const double dx = rod_length / (double)elements;
const double period = 60.0 / spm;
const double dt = period / (double)(outputs->point_count - 1);
if (!inputs->has_variable_rod || inputs->rod_node_count < nodes) {
inputs->has_variable_rod = 1;
inputs->rod_node_count = nodes;
const double E = 2.05e11;
const double rho = 7850.0;
const double taper = solver_input_or_default(inputs->taper_factor, 1.0);
const double area = solver_clamp((0.00042 + inputs->rod_friction_coefficient * 0.00008) * taper, 0.00025, 0.0009);
for (int i = 0; i < nodes; i++) {
inputs->area_m2[i] = area;
inputs->modulus_pa[i] = E;
inputs->density_kg_m3[i] = rho;
}
}
solver_trajectory_preprocess(inputs, elements, rod_length);
double kappa_mean = 0.0;
for (int i = 0; i < nodes; i++) {
kappa_mean += inputs->node_curvature[i];
}
kappa_mean /= (double)nodes;
const double mu_geo =
inputs->rod_friction_coefficient *
(inputs->molded_guide_mu_scale + inputs->wheeled_guide_mu_scale + inputs->other_guide_mu_scale) / 3.0;
double friction_geom = 1.0 + solver_clamp(mu_geo * kappa_mean * rod_length * 0.08, 0.0, 1.2);
if (!inputs->geometry_valid) {
friction_geom = solver_clamp(inputs->trajectory_friction_multiplier, 1.0, 1.8);
}
const double pump_friction = solver_clamp(inputs->pump_friction * friction_geom, 0.0, 7000.0);
const double stuffing_friction = solver_clamp(inputs->stuffing_box_friction * friction_geom, 0.0, 3500.0);
double side_nodes[SOLVER_MAX_NODES];
double fric_nodes[SOLVER_MAX_NODES];
memset(side_nodes, 0, sizeof(side_nodes));
memset(fric_nodes, 0, sizeof(fric_nodes));
const double k_pump = solver_clamp(1.8e5 + depth * 180.0, 1.2e5, 9.0e5);
const double c_pump = solver_clamp(220.0 + inputs->rod_friction_coefficient * 800.0, 120.0, 1200.0);
double baseline = depth * 2.05 + 3900.0;
for (int i = 0; i < elements; i++) {
const double rho_i = inputs->density_kg_m3[i];
const double A_i = 0.5 * (inputs->area_m2[i] + inputs->area_m2[i + 1]);
const double inc = 0.5 * (inputs->node_inc_rad[i] + inputs->node_inc_rad[i + 1]);
baseline += 0.5 * (rho_i - inputs->fluid_density_kg_m3) * A_i * inputs->gravity * cos(inc) * dx;
}
double u[MAX_FEA_NODES];
double v[MAX_FEA_NODES];
double a_vec[MAX_FEA_NODES];
memset(u, 0, sizeof(u));
memset(v, 0, sizeof(v));
memset(a_vec, 0, sizeof(a_vec));
const double surf_amp = solver_clamp(depth * 0.00045, 0.4, 2.6);
const double omega = 2.0 * M_PI / period;
solver_init_output_ranges(outputs);
if (inputs->workflow == 1 && inputs->surface_count >= 40) {
double pos_acc = 0.0;
double load_acc = 0.0;
for (int i = 0; i < inputs->surface_count; i++) {
pos_acc += inputs->surface_position_m[i];
load_acc += inputs->surface_load_n[i];
}
const double pm = pos_acc / (double)inputs->surface_count;
const double lm = load_acc / (double)inputs->surface_count;
const int steps = inputs->surface_count;
double dtd = dt;
if (inputs->surface_has_time && steps > 1) {
dtd = (inputs->surface_time_s[steps - 1] - inputs->surface_time_s[0]) / (double)(steps - 1);
}
for (int n = 0; n < steps; n++) {
const double u0 = inputs->surface_position_m[n] - pm;
const double v0 = (n + 1 < steps)
? (inputs->surface_position_m[n + 1] - inputs->surface_position_m[n]) / dtd
: (inputs->surface_position_m[0] - inputs->surface_position_m[steps - 1]) / dtd;
const double target_tension =
(inputs->surface_load_n[n] - lm) - baseline - stuffing_friction * solver_signum(v0);
double u_snap[MAX_FEA_NODES];
double v_snap[MAX_FEA_NODES];
double a_snap[MAX_FEA_NODES];
memcpy(u_snap, u, sizeof(u_snap));
memcpy(v_snap, v, sizeof(v_snap));
memcpy(a_snap, a_vec, sizeof(a_snap));
double lo = -1.2e6;
double hi = 1.2e6;
double fb = 0.0;
double pol_t = 0.0;
double dh_t = 0.0;
for (int it = 0; it < 24; it++) {
fb = 0.5 * (lo + hi);
memcpy(u, u_snap, sizeof(u_snap));
memcpy(v, v_snap, sizeof(v_snap));
memcpy(a_vec, a_snap, sizeof(a_snap));
if (fea_newmark_step(inputs, nodes, elements, dx, dtd, k_pump, c_pump, pump_friction, stuffing_friction, u0, v0, fb, u, v,
a_vec, &pol_t, &dh_t) != 0) {
return -3;
}
const double err = pol_t - target_tension;
if (fabs(err) < 40.0) {
break;
}
if (err > 0) {
hi = fb;
} else {
lo = fb;
}
}
const double polished = baseline + pol_t + stuffing_friction * solver_signum(v0);
const double downhole = baseline + dh_t;
outputs->position[n] = u0;
outputs->polished_load[n] = polished;
outputs->downhole_load[n] = downhole;
outputs->pump_position_m[n] = u[nodes - 1];
outputs->polished_stress_pa[n] = pol_t / fmax(0.5 * (inputs->area_m2[0] + inputs->area_m2[1]), 1e-12);
const double side_top = solver_compute_side_load_node(inputs, pol_t, 0, dx);
const double vloc = (n > 0) ? (outputs->pump_position_m[n] - outputs->pump_position_m[n - 1]) / dtd : v0;
outputs->side_load_profile_n[n] = side_top;
side_nodes[0] = side_top;
fric_nodes[0] = solver_compute_friction_node(inputs, side_top, vloc, 0);
side_nodes[nodes - 1] = solver_compute_side_load_node(inputs, dh_t, nodes - 1, dx);
fric_nodes[nodes - 1] = solver_compute_friction_node(inputs, side_nodes[nodes - 1], vloc, nodes - 1);
solver_update_output_ranges(outputs, polished, downhole);
solver_valve_state_step(inputs, outputs, n, outputs->pump_position_m[n], vloc, downhole);
}
outputs->point_count = steps;
for (int i = 0; i < steps; i++) {
if (i == 0) {
outputs->pump_velocity_m_s[i] = 0.0;
} else {
outputs->pump_velocity_m_s[i] = (outputs->pump_position_m[i] - outputs->pump_position_m[i - 1]) / dtd;
}
}
if (steps > 1) {
outputs->pump_velocity_m_s[0] = outputs->pump_velocity_m_s[1];
}
outputs->gas_interference =
(inputs->pump_fillage_option == 2 || (inputs->percent_pump_fillage > 0.0 && inputs->percent_pump_fillage < 88.0)) ? 1 : 0;
solver_fill_profiles(inputs, outputs, nodes, rod_length, side_nodes, fric_nodes);
if (inputs->enable_fourier_baseline) {
(void)solver_compute_fourier_baseline(inputs, outputs);
}
return 0;
}
for (int n = 0; n < outputs->point_count; n++) {
const double t = n * dt;
const double u0 = surf_amp * sin(omega * t);
const double v0 = surf_amp * omega * cos(omega * t);
double pol = 0.0;
double dhh = 0.0;
if (fea_newmark_step(inputs, nodes, elements, dx, dt, k_pump, c_pump, pump_friction, stuffing_friction, u0, v0, 0.0, u, v, a_vec,
&pol, &dhh) != 0) {
solver_add_warning(outputs, "FEA linear solve singular; aborting");
return -2;
}
const double polished = baseline + pol + stuffing_friction * solver_signum(v0);
const double downhole = baseline + dhh;
outputs->position[n] = u0;
outputs->polished_load[n] = polished;
outputs->downhole_load[n] = downhole;
outputs->pump_position_m[n] = u[nodes - 1];
outputs->polished_stress_pa[n] = pol / fmax(0.5 * (inputs->area_m2[0] + inputs->area_m2[1]), 1e-12);
const double side_top = solver_compute_side_load_node(inputs, pol, 0, dx);
outputs->side_load_profile_n[n] = side_top;
side_nodes[0] = side_top;
fric_nodes[0] = solver_compute_friction_node(inputs, side_top, v0, 0);
side_nodes[nodes - 1] = solver_compute_side_load_node(inputs, dhh, nodes - 1, dx);
fric_nodes[nodes - 1] = solver_compute_friction_node(inputs, side_nodes[nodes - 1], v0, nodes - 1);
solver_update_output_ranges(outputs, polished, downhole);
solver_valve_state_step(inputs, outputs, n, outputs->pump_position_m[n], v0, downhole);
}
for (int i = 0; i < outputs->point_count; i++) {
if (i == 0) {
outputs->pump_velocity_m_s[i] = 0.0;
} else {
outputs->pump_velocity_m_s[i] = (outputs->pump_position_m[i] - outputs->pump_position_m[i - 1]) / dt;
}
}
if (outputs->point_count > 1) {
outputs->pump_velocity_m_s[0] = outputs->pump_velocity_m_s[1];
}
outputs->gas_interference =
(inputs->pump_fillage_option == 2 || (inputs->percent_pump_fillage > 0.0 && inputs->percent_pump_fillage < 88.0)) ? 1 : 0;
solver_fill_profiles(inputs, outputs, nodes, rod_length, side_nodes, fric_nodes);
if (inputs->enable_fourier_baseline) {
(void)solver_compute_fourier_baseline(inputs, outputs);
}
return 0;
}

57
solver-c/src/solver_fourier.c Normal file
View File

@@ -0,0 +1,57 @@
#include "solver.h"
#include "solver_internal.h"
#include <math.h>
#include <string.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
static void reconstruct_series(const double *in, int n, int harmonics, double *out) {
if (n <= 0) return;
double a0 = 0.0;
for (int i = 0; i < n; i++) a0 += in[i];
a0 /= (double)n;
for (int i = 0; i < n; i++) out[i] = a0;
const int hmax = (harmonics > n / 2) ? (n / 2) : harmonics;
for (int k = 1; k <= hmax; k++) {
double ak = 0.0;
double bk = 0.0;
for (int i = 0; i < n; i++) {
const double th = 2.0 * M_PI * (double)k * (double)i / (double)n;
ak += in[i] * cos(th);
bk += in[i] * sin(th);
}
ak *= 2.0 / (double)n;
bk *= 2.0 / (double)n;
for (int i = 0; i < n; i++) {
const double th = 2.0 * M_PI * (double)k * (double)i / (double)n;
out[i] += ak * cos(th) + bk * sin(th);
}
}
}
int solver_compute_fourier_baseline(const SolverInputs *inputs, SolverOutputs *outputs) {
(void)inputs;
if (!outputs || outputs->point_count <= 0) return -1;
const int n = outputs->point_count;
const int harmonics = solver_clamp(inputs->fourier_harmonics, 1, SOLVER_MAX_FOURIER_HARMONICS);
reconstruct_series(outputs->polished_load, n, harmonics, outputs->fourier_polished_load);
reconstruct_series(outputs->downhole_load, n, harmonics, outputs->fourier_downhole_load);
double rss_p = 0.0;
double rss_d = 0.0;
for (int i = 0; i < n; i++) {
const double ep = outputs->polished_load[i] - outputs->fourier_polished_load[i];
const double ed = outputs->downhole_load[i] - outputs->fourier_downhole_load[i];
rss_p += ep * ep;
rss_d += ed * ed;
}
outputs->fourier_harmonics_used = harmonics;
outputs->fourier_residual_rms_polished = sqrt(rss_p / (double)n);
outputs->fourier_residual_rms_downhole = sqrt(rss_d / (double)n);
return 0;
}

114
solver-c/src/trajectory.c Normal file
View File

@@ -0,0 +1,114 @@
#include "solver.h"
#include <math.h>
#include <string.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
static double interp1(const double *x, const double *y, int n, double xq) {
if (n <= 0) {
return 0.0;
}
if (xq <= x[0]) {
return y[0];
}
if (xq >= x[n - 1]) {
return y[n - 1];
}
for (int i = 0; i < n - 1; i++) {
if (xq >= x[i] && xq <= x[i + 1]) {
const double t = (xq - x[i]) / (x[i + 1] - x[i] + 1e-30);
return y[i] + t * (y[i + 1] - y[i]);
}
}
return y[n - 1];
}
/* Curvature proxy between survey stations (SPE-173970 style dogleg), then map to measured depth along rod */
static void survey_curvature_kappa(const SolverInputs *in, double *kappa_station, int *n_k) {
const int n = in->survey_station_count;
*n_k = 0;
if (n < 3) {
return;
}
for (int i = 1; i < n; i++) {
const double ds = fmax(in->survey_md_m[i] - in->survey_md_m[i - 1], 1e-6);
const double dInc = in->survey_inc_rad[i] - in->survey_inc_rad[i - 1];
const double dAzi = in->survey_azi_rad[i] - in->survey_azi_rad[i - 1];
const double incMid = 0.5 * (in->survey_inc_rad[i] + in->survey_inc_rad[i - 1]);
const double kappa = sqrt(dInc * dInc + pow(sin(incMid) * dAzi, 2)) / ds;
kappa_station[i - 1] = kappa;
}
*n_k = n - 1;
}
static void survey_tangent_vectors(const SolverInputs *in, double tx[], double ty[], double tz[]) {
for (int i = 0; i < in->survey_station_count; i++) {
const double inc = in->survey_inc_rad[i];
const double azi = in->survey_azi_rad[i];
tx[i] = sin(inc) * cos(azi);
ty[i] = sin(inc) * sin(azi);
tz[i] = cos(inc);
}
}
void solver_trajectory_preprocess(SolverInputs *inputs, int nx, double rod_length_m) {
if (!inputs || nx <= 0 || rod_length_m <= 0.0) {
return;
}
const int nodes = nx + 1;
if (nodes > SOLVER_MAX_NODES) {
return;
}
memset(inputs->node_curvature, 0, sizeof(inputs->node_curvature));
memset(inputs->node_inc_rad, 0, sizeof(inputs->node_inc_rad));
memset(inputs->node_azi_rad, 0, sizeof(inputs->node_azi_rad));
memset(inputs->node_side_load_n, 0, sizeof(inputs->node_side_load_n));
if (inputs->survey_station_count < 3) {
inputs->geometry_valid = 0;
return;
}
double kappa_st[SOLVER_MAX_SURVEY];
double tx[SOLVER_MAX_SURVEY];
double ty[SOLVER_MAX_SURVEY];
double tz[SOLVER_MAX_SURVEY];
int nk = 0;
survey_curvature_kappa(inputs, kappa_st, &nk);
survey_tangent_vectors(inputs, tx, ty, tz);
if (nk <= 0) {
inputs->geometry_valid = 0;
return;
}
/* kappa_station[i] lives between md[i] and md[i+1]; map midpoint md */
double md_mid[SOLVER_MAX_SURVEY];
for (int i = 0; i < nk; i++) {
md_mid[i] = 0.5 * (inputs->survey_md_m[i] + inputs->survey_md_m[i + 1]);
}
const double md_top = inputs->survey_md_m[0];
const double md_bot = inputs->survey_md_m[inputs->survey_station_count - 1];
for (int j = 0; j < nodes; j++) {
const double s = rod_length_m * (double)j / (double)nx;
const double md_target = md_top + (md_bot - md_top) * (s / fmax(rod_length_m, 1e-6));
const double kappa = interp1(md_mid, kappa_st, nk, md_target);
const double tix = interp1(inputs->survey_md_m, tx, inputs->survey_station_count, md_target);
const double tiy = interp1(inputs->survey_md_m, ty, inputs->survey_station_count, md_target);
const double tiz = interp1(inputs->survey_md_m, tz, inputs->survey_station_count, md_target);
const double tnorm = fmax(sqrt(tix * tix + tiy * tiy + tiz * tiz), 1e-12);
const double nxv = tix / tnorm;
const double nyv = tiy / tnorm;
const double nzv = tiz / tnorm;
const double inc = acos(fmax(-1.0, fmin(1.0, nzv)));
const double azi = atan2(nyv, nxv);
inputs->node_curvature[j] = kappa;
inputs->node_inc_rad[j] = inc;
inputs->node_azi_rad[j] = azi;
}
inputs->geometry_valid = 1;
}

252
solver-c/tests/test_solver.c Normal file
View File

@@ -0,0 +1,252 @@
#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;
}
/* 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;
}
printf("solver-c tests passed\n");
return 0;
}