diff --git a/dev/C/Makefile b/dev/C/Makefile new file mode 100644 index 0000000..57a7a67 --- /dev/null +++ b/dev/C/Makefile @@ -0,0 +1,15 @@ +CC=gcc +PYTHON_VERSION=$(shell python3 -c "import sys; print(f'{sys.version_info.major}.{sys.version_info.minor}')") +PYTHON_INCLUDES=$(shell python3-config --includes) +PYTHON_LDFLAGS=$(shell python3-config --ldflags) +PYTHON_LIBS=$(shell python3-config --libs) + +# Add -lpython3.11 explicitly and move LDFLAGS to end of command +CFLAGS=-Wall -O3 $(PYTHON_INCLUDES) +LDFLAGS=$(PYTHON_LDFLAGS) $(PYTHON_LIBS) -lpython$(PYTHON_VERSION) + +fluids_test: fluids_test.c + $(CC) $(CFLAGS) $^ -o $@ $(LDFLAGS) + +clean: + rm -f fluids_test diff --git a/dev/C/fluids_test.c b/dev/C/fluids_test.c new file mode 100644 index 0000000..64c5b7f --- /dev/null +++ b/dev/C/fluids_test.c @@ -0,0 +1,430 @@ +#define PY_SSIZE_T_CLEAN +#include +#include +#include +#include +#include + +// Helper function to initialize Python and import fluids +PyObject* init_fluids(void) { + Py_Initialize(); + PyObject* module = PyImport_ImportModule("fluids"); + if (!module) { + PyErr_Print(); + fprintf(stderr, "Failed to import fluids module\n"); + return NULL; + } + return module; +} + +// Helper to measure time in microseconds +long long microseconds_now(void) { + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + return (long long)ts.tv_sec * 1000000LL + ts.tv_nsec / 1000LL; +} + +void test_atmosphere(PyObject* fluids) { + printf("\nTesting atmosphere at 5000m elevation:\n"); + + // Create ATMOSPHERE_1976 instance + PyObject* atm_class = PyObject_GetAttrString(fluids, "ATMOSPHERE_1976"); + if (!atm_class) { + PyErr_Print(); + return; + } + + PyObject* args = PyTuple_New(0); + PyObject* kwargs = PyDict_New(); + PyDict_SetItemString(kwargs, "Z", PyFloat_FromDouble(5000)); + + PyObject* atm = PyObject_Call(atm_class, args, kwargs); + if (!atm) { + PyErr_Print(); + Py_DECREF(args); + Py_DECREF(kwargs); + Py_DECREF(atm_class); + return; + } + + // Test various properties + PyObject* temp = PyObject_GetAttrString(atm, "T"); + PyObject* pressure = PyObject_GetAttrString(atm, "P"); + PyObject* density = PyObject_GetAttrString(atm, "rho"); + PyObject* gravity = PyObject_GetAttrString(atm, "g"); + PyObject* viscosity = PyObject_GetAttrString(atm, "mu"); + PyObject* conductivity = PyObject_GetAttrString(atm, "k"); + PyObject* sonic = PyObject_GetAttrString(atm, "v_sonic"); + + printf("✓ Temperature: %.4f\n", PyFloat_AsDouble(temp)); + printf("✓ Pressure: %.4f\n", PyFloat_AsDouble(pressure)); + printf("✓ Density: %.6f\n", PyFloat_AsDouble(density)); + printf("✓ Gravity: %.6f\n", PyFloat_AsDouble(gravity)); + printf("✓ Viscosity: %.6e\n", PyFloat_AsDouble(viscosity)); + printf("✓ Thermal conductivity: %.4f\n", PyFloat_AsDouble(conductivity)); + printf("✓ Sonic velocity: %.4f\n", PyFloat_AsDouble(sonic)); + + // Test static gravity method + PyObject* gravity_method = PyObject_GetAttrString(atm_class, "gravity"); + PyObject* gravity_args = PyTuple_New(0); + PyObject* gravity_kwargs = PyDict_New(); + PyDict_SetItemString(gravity_kwargs, "Z", PyFloat_FromDouble(1E5)); + + PyObject* high_gravity = PyObject_Call(gravity_method, gravity_args, gravity_kwargs); + printf("✓ High altitude gravity: %.6f\n", PyFloat_AsDouble(high_gravity)); + + // Cleanup + Py_DECREF(temp); + Py_DECREF(pressure); + Py_DECREF(density); + Py_DECREF(gravity); + Py_DECREF(viscosity); + Py_DECREF(conductivity); + Py_DECREF(sonic); + Py_DECREF(gravity_method); + Py_DECREF(gravity_args); + Py_DECREF(gravity_kwargs); + Py_DECREF(high_gravity); + Py_DECREF(atm); + Py_DECREF(args); + Py_DECREF(kwargs); + Py_DECREF(atm_class); +} + +void test_expanded_reynolds(PyObject* fluids) { + printf("\nTesting Reynolds number calculations:\n"); + + // Get Reynolds function + PyObject* reynolds_func = PyObject_GetAttrString(fluids, "Reynolds"); + if (!reynolds_func) { + PyErr_Print(); + return; + } + + // Test with density and viscosity + PyObject* args1 = PyTuple_New(0); + PyObject* kwargs1 = PyDict_New(); + PyDict_SetItemString(kwargs1, "V", PyFloat_FromDouble(2.5)); + PyDict_SetItemString(kwargs1, "D", PyFloat_FromDouble(0.25)); + PyDict_SetItemString(kwargs1, "rho", PyFloat_FromDouble(1.1613)); + PyDict_SetItemString(kwargs1, "mu", PyFloat_FromDouble(1.9E-5)); + + PyObject* re1 = PyObject_Call(reynolds_func, args1, kwargs1); + double re1_val = PyFloat_AsDouble(re1); + printf("✓ Re (with rho, mu): %.4f\n", re1_val); + assert(fabs(re1_val - 38200.6579) < 0.1); + + // Test with kinematic viscosity + PyObject* args2 = PyTuple_New(0); + PyObject* kwargs2 = PyDict_New(); + PyDict_SetItemString(kwargs2, "V", PyFloat_FromDouble(2.5)); + PyDict_SetItemString(kwargs2, "D", PyFloat_FromDouble(0.25)); + PyDict_SetItemString(kwargs2, "nu", PyFloat_FromDouble(1.636e-05)); + + PyObject* re2 = PyObject_Call(reynolds_func, args2, kwargs2); + double re2_val = PyFloat_AsDouble(re2); + printf("✓ Re (with nu): %.4f\n", re2_val); + assert(fabs(re2_val - 38202.934) < 0.1); + + // Cleanup + Py_DECREF(reynolds_func); + Py_DECREF(args1); + Py_DECREF(kwargs1); + Py_DECREF(re1); + Py_DECREF(args2); + Py_DECREF(kwargs2); + Py_DECREF(re2); +} + +void test_psd(PyObject* fluids) { + printf("\nTesting particle size distribution functionality:\n"); + + // Get particle_size_distribution module + PyObject* psd_module = PyObject_GetAttrString(fluids, "particle_size_distribution"); + if (!psd_module) { + PyErr_Print(); + return; + } + + // Create discrete PSD + PyObject* psd_class = PyObject_GetAttrString(psd_module, "ParticleSizeDistribution"); + if (!psd_class) { + PyErr_Print(); + Py_DECREF(psd_module); + return; + } + + // Create lists for ds and numbers + PyObject* ds_list = PyList_New(15); + double ds[] = {240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, + 2141, 2676, 3345, 4181, 5226, 6532}; + for (int i = 0; i < 15; i++) { + PyList_SetItem(ds_list, i, PyFloat_FromDouble(ds[i])); + } + + PyObject* numbers_list = PyList_New(14); + double numbers[] = {65, 119, 232, 410, 629, 849, 990, 981, 825, + 579, 297, 111, 21, 1}; + for (int i = 0; i < 14; i++) { + PyList_SetItem(numbers_list, i, PyFloat_FromDouble(numbers[i])); + } + + PyObject* args = PyTuple_New(0); + PyObject* kwargs = PyDict_New(); + PyDict_SetItemString(kwargs, "ds", ds_list); + PyDict_SetItemString(kwargs, "fractions", numbers_list); + PyDict_SetItemString(kwargs, "order", PyLong_FromLong(0)); + + PyObject* psd = PyObject_Call(psd_class, args, kwargs); + printf("✓ Created discrete PSD\n"); + + // Test mean sizes + PyObject* mean_size_method = PyObject_GetAttrString(psd, "mean_size"); + PyObject* mean_args = PyTuple_Pack(2, PyLong_FromLong(2), PyLong_FromLong(1)); + PyObject* d21 = PyObject_Call(mean_size_method, mean_args, NULL); + printf("✓ Size-weighted mean diameter: %.4f\n", PyFloat_AsDouble(d21)); + + PyObject* mean_args2 = PyTuple_Pack(2, PyLong_FromLong(1), PyLong_FromLong(0)); + PyObject* d10 = PyObject_Call(mean_size_method, mean_args2, NULL); + printf("✓ Arithmetic mean diameter: %.4f\n", PyFloat_AsDouble(d10)); + + // Test percentile calculations + PyObject* dn_method = PyObject_GetAttrString(psd, "dn"); + PyObject* d10_args = PyTuple_Pack(1, PyFloat_FromDouble(0.1)); + PyObject* d90_args = PyTuple_Pack(1, PyFloat_FromDouble(0.9)); + + PyObject* d10_percentile = PyObject_Call(dn_method, d10_args, NULL); + PyObject* d90_percentile = PyObject_Call(dn_method, d90_args, NULL); + + printf("✓ D10: %.4f\n", PyFloat_AsDouble(d10_percentile)); + printf("✓ D90: %.4f\n", PyFloat_AsDouble(d90_percentile)); + + // Cleanup + Py_DECREF(psd_module); + Py_DECREF(psd_class); + Py_DECREF(ds_list); + Py_DECREF(numbers_list); + Py_DECREF(args); + Py_DECREF(kwargs); + Py_DECREF(psd); + Py_DECREF(mean_size_method); + Py_DECREF(mean_args); + Py_DECREF(d21); + Py_DECREF(mean_args2); + Py_DECREF(d10); + Py_DECREF(dn_method); + Py_DECREF(d10_args); + Py_DECREF(d90_args); + Py_DECREF(d10_percentile); + Py_DECREF(d90_percentile); +} + +void test_fluids(PyObject* fluids) { + printf("Running fluids tests from C...\n"); + + // Get version + PyObject* version = PyObject_GetAttrString(fluids, "__version__"); + if (version) { + const char* version_str = PyUnicode_AsUTF8(version); + printf("✓ Successfully imported fluids\n"); + printf("✓ Fluids version: %s\n", version_str); + Py_DECREF(version); + } + + // Test Reynolds number calculation + PyObject* reynolds_func = PyObject_GetAttrString(fluids, "Reynolds"); + if (reynolds_func) { + PyObject* args = PyTuple_New(0); + PyObject* kwargs = PyDict_New(); + PyDict_SetItemString(kwargs, "V", PyFloat_FromDouble(2.5)); + PyDict_SetItemString(kwargs, "D", PyFloat_FromDouble(0.1)); + PyDict_SetItemString(kwargs, "rho", PyFloat_FromDouble(1000)); + PyDict_SetItemString(kwargs, "mu", PyFloat_FromDouble(0.001)); + + PyObject* result = PyObject_Call(reynolds_func, args, kwargs); + if (result) { + double re = PyFloat_AsDouble(result); + printf("✓ Reynolds number calculation successful: %f\n", re); + assert(re > 0); + Py_DECREF(result); + } + + Py_DECREF(args); + Py_DECREF(kwargs); + Py_DECREF(reynolds_func); + } +} + +void benchmark_fluids(PyObject* fluids) { + printf("\nRunning benchmarks:\n"); + + // Get friction_factor function + PyObject* friction_func = PyObject_GetAttrString(fluids, "friction_factor"); + if (!friction_func) { + PyErr_Print(); + return; + } + + // Get TANK class + PyObject* tank_class = PyObject_GetAttrString(fluids, "TANK"); + if (!tank_class) { + PyErr_Print(); + Py_DECREF(friction_func); + return; + } + + // Benchmark friction_factor + printf("\nBenchmarking friction_factor:\n"); + long long start = microseconds_now(); + + for(int i = 0; i < 1000000; i++) { + PyObject* args = PyTuple_New(0); + PyObject* kwargs = PyDict_New(); + PyDict_SetItemString(kwargs, "Re", PyFloat_FromDouble(1e5)); + PyDict_SetItemString(kwargs, "eD", PyFloat_FromDouble(0.0001)); + + PyObject* result = PyObject_Call(friction_func, args, kwargs); + + Py_DECREF(result); + Py_DECREF(args); + Py_DECREF(kwargs); + } + + long long duration = microseconds_now() - start; + printf("Time for 1e6 friction_factor calls: %lld microseconds\n", duration); + printf("Average time per call: %.6f microseconds\n", duration / 1000000.0); + + // Benchmark TANK creation + printf("\nBenchmarking TANK creation:\n"); + start = microseconds_now(); + + for(int i = 0; i < 1000; i++) { + PyObject* args = PyTuple_New(0); + PyObject* kwargs = PyDict_New(); + PyDict_SetItemString(kwargs, "L", PyFloat_FromDouble(3)); + PyDict_SetItemString(kwargs, "D", PyFloat_FromDouble(5)); + PyDict_SetItemString(kwargs, "horizontal", Py_False); + PyDict_SetItemString(kwargs, "sideA", PyUnicode_FromString("torispherical")); + PyDict_SetItemString(kwargs, "sideB", PyUnicode_FromString("torispherical")); + PyDict_SetItemString(kwargs, "sideA_f", PyFloat_FromDouble(1)); + PyDict_SetItemString(kwargs, "sideA_k", PyFloat_FromDouble(0.1)); + PyDict_SetItemString(kwargs, "sideB_f", PyFloat_FromDouble(1)); + PyDict_SetItemString(kwargs, "sideB_k", PyFloat_FromDouble(0.1)); + + PyObject* tank = PyObject_Call(tank_class, args, kwargs); + + Py_DECREF(tank); + Py_DECREF(args); + Py_DECREF(kwargs); + } + + duration = microseconds_now() - start; + printf("Average time per tank creation: %.6f microseconds\n", duration / 1000.0); + + Py_DECREF(friction_func); + Py_DECREF(tank_class); +} + +void test_tank(PyObject* fluids) { + printf("\nTesting tank calculations:\n"); + + // Get TANK class + PyObject* tank_class = PyObject_GetAttrString(fluids, "TANK"); + if (!tank_class) { + PyErr_Print(); + return; + } + + // Test basic tank creation + PyObject* args1 = PyTuple_New(0); + PyObject* kwargs1 = PyDict_New(); + PyDict_SetItemString(kwargs1, "V", PyFloat_FromDouble(10)); + PyDict_SetItemString(kwargs1, "L_over_D", PyFloat_FromDouble(0.7)); + PyDict_SetItemString(kwargs1, "sideB", PyUnicode_FromString("conical")); + PyDict_SetItemString(kwargs1, "horizontal", Py_False); + + PyObject* T1 = PyObject_Call(tank_class, args1, kwargs1); + if (!T1) { + PyErr_Print(); + goto cleanup1; + } + + PyObject* length = PyObject_GetAttrString(T1, "L"); + PyObject* diameter = PyObject_GetAttrString(T1, "D"); + printf("✓ Tank length: %.6f\n", PyFloat_AsDouble(length)); + printf("✓ Tank diameter: %.6f\n", PyFloat_AsDouble(diameter)); + + // Test torispherical tank + PyObject* args2 = PyTuple_New(0); + PyObject* kwargs2 = PyDict_New(); + PyDict_SetItemString(kwargs2, "L", PyFloat_FromDouble(3)); + PyDict_SetItemString(kwargs2, "D", PyFloat_FromDouble(5)); + PyDict_SetItemString(kwargs2, "horizontal", Py_False); + PyDict_SetItemString(kwargs2, "sideA", PyUnicode_FromString("torispherical")); + PyDict_SetItemString(kwargs2, "sideB", PyUnicode_FromString("torispherical")); + PyDict_SetItemString(kwargs2, "sideA_f", PyFloat_FromDouble(1)); + PyDict_SetItemString(kwargs2, "sideA_k", PyFloat_FromDouble(0.1)); + PyDict_SetItemString(kwargs2, "sideB_f", PyFloat_FromDouble(1)); + PyDict_SetItemString(kwargs2, "sideB_k", PyFloat_FromDouble(0.1)); + + PyObject* DIN = PyObject_Call(tank_class, args2, kwargs2); + if (!DIN) { + PyErr_Print(); + goto cleanup2; + } + + // Get max height + PyObject* h_max = PyObject_GetAttrString(DIN, "h_max"); + printf("✓ Tank max height: %.6f\n", PyFloat_AsDouble(h_max)); + + // Test h_from_V method + PyObject* h_from_V = PyObject_GetAttrString(DIN, "h_from_V"); + PyObject* h_args = PyTuple_Pack(1, PyFloat_FromDouble(40)); + PyObject* h_result = PyObject_CallObject(h_from_V, h_args); + printf("✓ Height at V=40: %.6f\n", PyFloat_AsDouble(h_result)); + + // Test V_from_h method + PyObject* V_from_h = PyObject_GetAttrString(DIN, "V_from_h"); + PyObject* v_args = PyTuple_Pack(1, PyFloat_FromDouble(4.1)); + PyObject* v_result = PyObject_CallObject(V_from_h, v_args); + printf("✓ Volume at h=4.1: %.5f\n", PyFloat_AsDouble(v_result)); + + // Cleanup + Py_DECREF(v_result); + Py_DECREF(v_args); + Py_DECREF(V_from_h); + Py_DECREF(h_result); + Py_DECREF(h_args); + Py_DECREF(h_from_V); + Py_DECREF(h_max); + Py_DECREF(DIN); +cleanup2: + Py_DECREF(args2); + Py_DECREF(kwargs2); + Py_DECREF(diameter); + Py_DECREF(length); + Py_DECREF(T1); +cleanup1: + Py_DECREF(args1); + Py_DECREF(kwargs1); + Py_DECREF(tank_class); +} + +int main() { + PyObject* fluids = init_fluids(); + if (!fluids) { + return 1; + } + + test_fluids(fluids); + test_atmosphere(fluids); + test_tank(fluids); + test_expanded_reynolds(fluids); + test_psd(fluids); + benchmark_fluids(fluids); + + Py_DECREF(fluids); + Py_Finalize(); + printf("\nAll tests completed!\n"); + return 0; +} diff --git a/dev/D/dub.json b/dev/D/dub.json new file mode 100644 index 0000000..53a9a8c --- /dev/null +++ b/dev/D/dub.json @@ -0,0 +1,16 @@ +{ + "name": "pyd-test", + "description": "A test PyD project", + "dependencies": { + "pyd": "~>0.14.4" + }, + "targetType": "executable", + "configurations": [ + { + "name": "python311", + "subConfigurations": { + "pyd": "python311" + } + } + ] +} diff --git a/dev/D/dub.selections.json b/dev/D/dub.selections.json new file mode 100644 index 0000000..ea88456 --- /dev/null +++ b/dev/D/dub.selections.json @@ -0,0 +1,6 @@ +{ + "fileVersion": 1, + "versions": { + "pyd": "0.14.4" + } +} diff --git a/dev/D/source/app.d b/dev/D/source/app.d new file mode 100644 index 0000000..2d56d18 --- /dev/null +++ b/dev/D/source/app.d @@ -0,0 +1,304 @@ +// import std.stdio; +// import pyd.pyd; +// import pyd.embedded; +// import std.format; +// +// // Initialize Python before using it +// shared static this() { +// py_init(); +// } +// +// // Helper to run Python expressions and handle errors +// auto pyEval(T)(string expr, string namespace = "") { +// try { +// return py_eval!T(expr, namespace); +// } catch (Exception e) { +// writeln("Error evaluating: ", expr); +// writeln("Error: ", e.msg); +// throw e; +// } +// } +// +// void testFluids() { +// writeln("\nTesting basic fluids functionality..."); +// +// // Import fluids and create a PydObject reference +// PydObject fluidsModule = py_import("fluids"); +// writeln("✓ Successfully imported fluids"); +// +// // Get version using attribute access +// string fluidsVersion = fluidsModule.__version__.to_d!string(); +// writeln("✓ Fluids version: ", fluidsVersion); +// +// double Re = pyEval!double("Reynolds(V=2.5, D=0.1, rho=1000, mu=0.001)", "fluids"); +// writeln("✓ Reynolds number calculation: ", Re); +// assert(Re > 0); +// +// // Test friction factor +// double fd = pyEval!double("friction_factor(Re=1e5, eD=0.0001)", "fluids"); +// writeln("✓ Friction factor: ", fd); +// assert(0 < fd && fd < 1); +// } +// +// void testAtmosphere() { +// writeln("\nTesting atmosphere calculations..."); +// +// PydObject fluidsModule = py_import("fluids"); +// PydObject atm = py_eval("ATMOSPHERE_1976(5000.0)", "fluids"); +// +// +// // Get properties +// double temp = atm.T.to_d!double(); +// double pressure = atm.P.to_d!double(); +// double density = atm.rho.to_d!double(); +// +// writeln("At 5000m elevation:"); +// writeln("✓ Temperature: ", temp); +// writeln("✓ Pressure: ", pressure); +// writeln("✓ Density: ", density); +// } +// +// void main() { +// try { +// writeln("Running fluids tests from D..."); +// testFluids(); +// testAtmosphere(); +// writeln("\nAll tests completed!"); +// } catch (Exception e) { +// writeln("Test failed with error: ", e.msg); +// } +// } + + + + +import std.stdio; +import pyd.pyd; +import pyd.embedded; +import std.format; +import std.math; +import std.array; +import core.time; +import std.datetime.stopwatch; + +// Initialize Python before using it +shared static this() { + py_init(); +} + +// Helper to run Python expressions and handle errors +auto pyEval(T)(string expr, string namespace = "") { + try { + return py_eval!T(expr, namespace); + } catch (Exception e) { + writeln("Error evaluating: ", expr); + writeln("Error: ", e.msg); + throw e; + } +} + +void testFluids() { + writeln("\nTesting basic fluids functionality..."); + + // Import fluids and create a PydObject reference + PydObject fluidsModule = py_import("fluids"); + writeln("✓ Successfully imported fluids"); + + // Get version using attribute access + string fluidsVersion = fluidsModule.__version__.to_d!string(); + writeln("✓ Fluids version: ", fluidsVersion); + + double Re = pyEval!double("Reynolds(V=2.5, D=0.1, rho=1000, mu=0.001)", "fluids"); + writeln("✓ Reynolds number calculation: ", Re); + assert(Re > 0); + + // Test friction factor + double fd = pyEval!double("friction_factor(Re=1e5, eD=0.0001)", "fluids"); + writeln("✓ Friction factor: ", fd); + assert(0 < fd && fd < 1); +} + +void testAtmosphere() { + writeln("\nTesting atmosphere calculations..."); + PydObject fluidsModule = py_import("fluids"); + PydObject atm = py_eval("ATMOSPHERE_1976(5000.0)", "fluids"); + + // Get properties + double temp = atm.T.to_d!double(); + double pressure = atm.P.to_d!double(); + double density = atm.rho.to_d!double(); + double gravity = atm.g.to_d!double(); + double viscosity = atm.mu.to_d!double(); + double thermalConductivity = atm.k.to_d!double(); + double sonicVelocity = atm.v_sonic.to_d!double(); + + writeln("At 5000m elevation:"); + writefln("✓ Temperature: %.4f", temp); + writefln("✓ Pressure: %.4f", pressure); + writefln("✓ Density: %.6f", density); + writefln("✓ Gravity: %.6f", gravity); + writefln("✓ Viscosity: %.6e", viscosity); + writefln("✓ Thermal conductivity: %.6f", thermalConductivity); + writefln("✓ Sonic velocity: %.4f", sonicVelocity); + + // Test static methods + double gHigh = pyEval!double("ATMOSPHERE_1976.gravity(1E5)", "fluids"); + writefln("✓ High altitude gravity: %.6f", gHigh); + + double vSonic = pyEval!double("ATMOSPHERE_1976.sonic_velocity(300)", "fluids"); + writefln("✓ Sonic velocity at 300K: %.4f", vSonic); + + double mu400 = pyEval!double("ATMOSPHERE_1976.viscosity(400)", "fluids"); + writefln("✓ Viscosity at 400K: %.6e", mu400); + + double k400 = pyEval!double("ATMOSPHERE_1976.thermal_conductivity(400)", "fluids"); + writefln("✓ Thermal conductivity at 400K: %.6f", k400); +} + +void testTank() { + writeln("\nTesting tank calculations..."); + PydObject fluidsModule = py_import("fluids"); + + // Test basic tank creation + PydObject T1 = py_eval("TANK(V=10, L_over_D=0.7, sideB='conical', horizontal=False)", "fluids"); + writefln("✓ Tank length: %.6f", T1.L.to_d!double()); + writefln("✓ Tank diameter: %.6f", T1.D.to_d!double()); + + // Test ellipsoidal tank + PydObject tankEllip = py_eval( + "TANK(D=10, V=500, horizontal=False, sideA='ellipsoidal', sideB='ellipsoidal', sideA_a=1, sideB_a=1)", + "fluids" + ); + writefln("✓ Ellipsoidal tank L: %.6f", tankEllip.L.to_d!double()); + + // Test torispherical tank + PydObject DIN = py_eval( + "TANK(L=3, D=5, horizontal=False, sideA='torispherical', sideB='torispherical', " ~ + "sideA_f=1, sideA_k=0.1, sideB_f=1, sideB_k=0.1)", + "fluids" + ); + + writeln("✓ Tank representation: ", DIN.toString()); + writefln("✓ Height at V=40: %.6f", DIN.method("h_from_V", 40.0).to_d!double()); + writefln("✓ Volume at h=4.1: %.5f", DIN.method("V_from_h", 4.1).to_d!double()); + writefln("✓ Surface area at h=2.1: %.5f", DIN.method("SA_from_h", 2.1).to_d!double()); + +} + +void testReynolds() { + writeln("\nTesting Reynolds number calculations:"); + + // Test with density and viscosity + double Re1 = pyEval!double("Reynolds(V=2.5, D=0.25, rho=1.1613, mu=1.9E-5)", "fluids"); + writefln("✓ Re (with rho, mu): %.4f", Re1); + assert(abs(Re1 - 38200.6579) < 0.1); + + // Test with kinematic viscosity + double Re2 = pyEval!double("Reynolds(V=2.5, D=0.25, nu=1.636e-05)", "fluids"); + writefln("✓ Re (with nu): %.4f", Re2); + assert(abs(Re2 - 38202.934) < 0.1); +} + +void testPSD() { + writeln("\nTesting particle size distribution functionality:"); + + // Create arrays for discrete PSD + double[] ds = [240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532]; + double[] numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1]; + + // Create Python lists from D arrays + string dsStr = format("[%(%s,%)]", ds); + string numbersStr = format("[%(%s,%)]", numbers); + + // Create the PSD object + PydObject particleDist = py_eval("particle_size_distribution", "fluids"); + PydObject psd = particleDist.ParticleSizeDistribution(ds, numbers, 0); + writeln("✓ Created discrete PSD"); + + // Test mean sizes + double d21 = psd.method("mean_size", 2, 1).to_d!double(); + writefln("✓ Size-weighted mean diameter: %.4f", d21); + + double d10 = psd.method("mean_size", 1, 0).to_d!double(); + writefln("✓ Arithmetic mean diameter: %.4f", d10); + + // Test percentile calculations + double d10Percentile = psd.method("dn", 0.1).to_d!double(); + double d90Percentile = psd.method("dn", 0.9).to_d!double(); + writefln("✓ D10: %.4f", d10Percentile); + writefln("✓ D90: %.4f", d90Percentile); + + // Test probability functions + double pdfVal = psd.method("pdf", 1000.0).to_d!double(); + double cdfVal = psd.method("cdf", 5000.0).to_d!double(); + writefln("✓ PDF at 1000: %.4e", pdfVal); + writefln("✓ CDF at 5000: %.6f", cdfVal); + + // Test lognormal distribution + PydObject psdLog = particleDist.PSDLognormal(0.5, 5e-6); + writeln("✓ Created lognormal PSD"); + + double vssa = psdLog.vssa.to_d!double(); + writefln("✓ Volume specific surface area: %.2f", vssa); + + // Calculate span using individual method calls + double dn90 = psdLog.method("dn", 0.9).to_d!double(); + double dn10 = psdLog.method("dn", 0.1).to_d!double(); + double span = dn90 - dn10; + writefln("✓ Span: %.4e", span); + + // Calculate ratio using individual method calls + double dn75 = psdLog.method("dn", 0.75).to_d!double(); + double dn25 = psdLog.method("dn", 0.25).to_d!double(); + double ratio7525 = dn75 / dn25; + writefln("✓ D75/D25 ratio: %.6f", ratio7525); +} + +void benchmarkFluids() { + writeln("\nRunning benchmarks:"); + + // Benchmark friction factor calculation + writeln("\nBenchmarking friction_factor:"); + auto sw = StopWatch(AutoStart.yes); + + foreach (i; 0..10000) { + pyEval!double("friction_factor(Re=1e5, eD=0.0001)", "fluids"); + } + + sw.stop(); + double elapsed = sw.peek().total!"usecs" / 1_000_000.0; + writefln("Time for 1e4 friction_factor calls: %.6f seconds", elapsed); + writefln("Average time per call: %.6f seconds", elapsed/10000); + + // Benchmark tank creation + writeln("\nBenchmarking TANK creation:"); + sw.reset(); + sw.start(); + + foreach (i; 0..1_000) { + py_eval( + "TANK(L=3, D=5, horizontal=False, sideA='torispherical', sideB='torispherical', " ~ + "sideA_f=1, sideA_k=0.1, sideB_f=1, sideB_k=0.1)", + "fluids" + ); + } + + sw.stop(); + elapsed = sw.peek().total!"usecs" / 1_000_000.0; + writefln("Average time per creation: %.6f seconds", elapsed/1_000); +} + +void main() { + try { + writeln("Running fluids tests from D..."); + testFluids(); + testAtmosphere(); + testTank(); + testReynolds(); + testPSD(); + benchmarkFluids(); + writeln("\nAll tests completed!"); + } catch (Exception e) { + writeln("Test failed with error: ", e.msg); + } +} diff --git a/dev/check_fluids_can_be_called_from_lisp.lisp b/dev/check_fluids_can_be_called_from_lisp.lisp new file mode 100755 index 0000000..4f0c5c7 --- /dev/null +++ b/dev/check_fluids_can_be_called_from_lisp.lisp @@ -0,0 +1,245 @@ +#!/usr/bin/sbcl --script + +;; Load quicklisp +#-quicklisp +(let ((quicklisp-init (merge-pathnames "quicklisp/setup.lisp" + (user-homedir-pathname)))) + (when (probe-file quicklisp-init) + (load quicklisp-init))) + +;; Load py4cl and ensure we're using Python 3 +(ql:quickload :py4cl :silent t) +(setf py4cl:*python-command* "python3") + +;; Import required modules +;; (py4cl:import-module "fluids" :as "fl") +(py4cl:import-module "fluids" :as "fluids") + +;; Check Python environment +(py4cl:python-exec " +import sys +print('Python path:', sys.executable) +try: + import numpy + print('NumPy found at:', numpy.__file__) + import fluids + print('fluids found at:', fluids.__file__) +except ImportError as e: + print('Error importing NumPy:', e) +") + + + +(defun test-fluids () + (handler-case + (progn + (format t "✓ Successfully imported fluids~%") + (format t "✓ Fluids version: ~A~%" (py4cl:python-eval "fluids.__version__")) + ;; Test Reynolds number calculation + (let ((re (py4cl:python-eval + (format nil "fluids.Reynolds(**{'V': ~F, 'D': ~F, 'rho': ~F, 'mu': ~F})" + 2.5 0.1 1000 0.001)))) + (format t "✓ Reynolds number calculation successful: ~A~%" re) + (assert (> re 0))) + + ;; Test friction factor calculation + (let ((fd (py4cl:python-eval + (format nil "fluids.friction_factor(**{'Re': ~E, 'eD': ~F})" + 1e5 0.0001)))) + (format t "✓ Friction factor calculation successful: ~A~%" fd) + (assert (and (> fd 0) (< fd 1)))) + + (format t "~%All basic tests completed successfully!~%")) + (error (e) + (format t "Error occurred: ~A~%" e) + (error e)))) + + + +(defun test-atmosphere () + (handler-case + (progn + ;; Create and store the atmosphere object in Python's namespace + (py4cl:python-exec + (format nil "atm = fluids.ATMOSPHERE_1976(**{'Z': ~F})" 5000)) + + (format t "~%Testing atmosphere at 5000m elevation:~%") + ;; Access properties using attributes + (format t "✓ Temperature: ~,4F~%" + (py4cl:python-eval "atm.T")) + (format t "✓ Pressure: ~,4F~%" + (py4cl:python-eval "atm.P")) + (format t "✓ Density: ~,6F~%" + (py4cl:python-eval "atm.rho")) + + ;; Test derived properties + (format t "✓ Gravity: ~,6F~%" + (py4cl:python-eval "atm.g")) + (format t "✓ Viscosity: ~,6E~%" + (py4cl:python-eval "atm.mu")) + (format t "✓ Thermal conductivity: ~,6F~%" + (py4cl:python-eval "atm.k")) + (format t "✓ Sonic velocity: ~,4F~%" + (py4cl:python-eval "atm.v_sonic")) + + ;; Test static methods + (let ((g-high (py4cl:python-eval + (format nil "fluids.ATMOSPHERE_1976.gravity(**{'Z': ~E})" 1e5)))) + (format t "✓ High altitude gravity: ~,6F~%" g-high)) + + (let ((v-sonic (py4cl:python-eval + (format nil "fluids.ATMOSPHERE_1976.sonic_velocity(**{'T': ~F})" 300)))) + (format t "✓ Sonic velocity at 300K: ~,4F~%" v-sonic)) + + (let ((mu-400 (py4cl:python-eval + (format nil "fluids.ATMOSPHERE_1976.viscosity(**{'T': ~F})" 400)))) + (format t "✓ Viscosity at 400K: ~,6E~%" mu-400)) + + (let ((k-400 (py4cl:python-eval + (format nil "fluids.ATMOSPHERE_1976.thermal_conductivity(**{'T': ~F})" 400)))) + (format t "✓ Thermal conductivity at 400K: ~,6F~%" k-400))) + (error (e) + (format t "Error in atmosphere tests: ~A~%" e) + (error e)))) + + +(defun test-tank () + (handler-case + (progn + ;; Test basic tank creation + (py4cl:python-exec + (format nil "t1 = fluids.TANK(**{'V': ~F, 'L_over_D': ~F, 'sideB': 'conical', 'horizontal': False})" + 10 0.7)) + (format t "~%Testing tank calculations:~%") + (format t "✓ Tank length: ~,6F~%" + (py4cl:python-eval "t1.L")) + (format t "✓ Tank diameter: ~,6F~%" + (py4cl:python-eval "t1.D")) + + ;; Test ellipsoidal tank + (py4cl:python-exec + (format nil "tank_ellip = fluids.TANK(**{'D': ~F, 'V': ~F, 'horizontal': False, + 'sideA': 'ellipsoidal', 'sideB': 'ellipsoidal', + 'sideA_a': ~F, 'sideB_a': ~F})" + 10 500 1 1)) + (format t "✓ Ellipsoidal tank L: ~,6F~%" + (py4cl:python-eval "tank_ellip.L")) + + ;; Test torispherical tank + (py4cl:python-exec + (format nil "din = fluids.TANK(**{'L': ~F, 'D': ~F, 'horizontal': False, + 'sideA': 'torispherical', 'sideB': 'torispherical', + 'sideA_f': ~F, 'sideA_k': ~F, 'sideB_f': ~F, 'sideB_k': ~F})" + 3 5 1 0.1 1 0.1)) + + (format t "✓ Tank representation: ~A~%" + (py4cl:python-eval "str(din)")) + (format t "✓ Tank max height: ~,6F~%" + (py4cl:python-eval "din.h_max")) + (format t "✓ Height at V=40: ~,6F~%" + (py4cl:python-eval "din.h_from_V(40)")) + (format t "✓ Volume at h=4.1: ~,5F~%" + (py4cl:python-eval "din.V_from_h(4.1)")) + (format t "✓ Surface area at h=2.1: ~,5F~%" + (py4cl:python-eval "din.SA_from_h(2.1)"))) + (error (e) + (format t "Error in tank tests: ~A~%" e) + (error e)))) + + +(defun test-psd () + (handler-case + (progn + (format t "~%Testing particle size distribution functionality:~%") + + ;; Create arrays in Python's namespace + (py4cl:python-exec " +ds = [240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532] +numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1]") + + ;; Create discrete PSD + (py4cl:python-exec " +psd = fluids.particle_size_distribution.ParticleSizeDistribution( + ds=ds, + fractions=numbers, + order=0 +)") + (format t "✓ Created discrete PSD~%") + + ;; Test mean sizes + (let ((d21 (py4cl:python-eval "psd.mean_size(2, 1)"))) + (format t "✓ Size-weighted mean diameter: ~,4F~%" d21) + (assert (< (abs (- d21 1857.788)) 0.1))) + + (let ((d10 (py4cl:python-eval "psd.mean_size(1, 0)"))) + (format t "✓ Arithmetic mean diameter: ~,4F~%" d10) + (assert (< (abs (- d10 1459.372)) 0.1))) + + ;; Test percentile calculations + (let ((d10-percentile (py4cl:python-eval "psd.dn(0.1)")) + (d90-percentile (py4cl:python-eval "psd.dn(0.9)"))) + (format t "✓ D10: ~,4F~%" d10-percentile) + (format t "✓ D90: ~,4F~%" d90-percentile)) + + ;; Test probability functions + (let ((pdf-val (py4cl:python-eval "psd.pdf(1000)")) + (cdf-val (py4cl:python-eval "psd.cdf(5000)"))) + (format t "✓ PDF at 1000: ~,4E~%" pdf-val) + (format t "✓ CDF at 5000: ~,6F~%" cdf-val)) + + ;; Test lognormal distribution + (py4cl:python-exec + (format nil "psd_log = fluids.particle_size_distribution.PSDLognormal(**{'s': ~F, 'd_characteristic': ~E})" + 0.5 5e-6)) + (format t "✓ Created lognormal PSD~%") + + (let ((vssa (py4cl:python-eval "psd_log.vssa"))) + (format t "✓ Volume specific surface area: ~,2F~%" vssa)) + + (let ((span (py4cl:python-eval "psd_log.dn(0.9) - psd_log.dn(0.1)"))) + (format t "✓ Span: ~,4E~%" span)) + + (let ((ratio-7525 (py4cl:python-eval "psd_log.dn(0.75)/psd_log.dn(0.25)"))) + (format t "✓ D75/D25 ratio: ~,6F~%" ratio-7525))) + (error (e) + (format t "Error in PSD tests: ~A~%" e) + (error e)))) + +(defun benchmark-fluids () + (format t "~%Running benchmarks:~%") + + ;; Benchmark friction factor calculation + (format t "~%Benchmarking friction_factor:~%") + (let ((t1 (get-internal-real-time))) + (dotimes (i 10000) + (py4cl:python-eval + (format nil "fluids.friction_factor(**{'Re': ~E, 'eD': ~F})" 1e5 0.0001))) + (let ((elapsed (/ (- (get-internal-real-time) t1) internal-time-units-per-second))) + (format t "Time for 1e6 friction_factor calls: ~,6F seconds~%" elapsed) + (format t "Average time per call: ~,6F seconds~%" (/ elapsed 10000)))) + + ;; Benchmark tank creation + (format t "~%Benchmarking TANK creation:~%") + (let ((t2 (get-internal-real-time))) + (dotimes (i 10000) + (py4cl:python-eval + (format nil "fluids.TANK(**{'L': ~F, 'D': ~F, 'horizontal': False, + 'sideA': 'torispherical', 'sideB': 'torispherical', + 'sideA_f': ~F, 'sideA_k': ~F, + 'sideB_f': ~F, 'sideB_k': ~F})" + 3 5 1 0.1 1 0.1))) + (let ((elapsed (/ (- (get-internal-real-time) t2) internal-time-units-per-second))) + (format t "Average time per creation: ~,6F seconds~%" (/ elapsed 10000))))) + + + ;; Run all tests +(format t "Running fluids tests from Common Lisp...~%") +(test-fluids) +(test-atmosphere) +(test-tank) +(test-psd) +(benchmark-fluids) +(format t "~%All tests completed!~%") + +;; Clean up +(py4cl:python-stop) diff --git a/dev/check_fluids_can_be_called_from_perl.pl b/dev/check_fluids_can_be_called_from_perl.pl new file mode 100755 index 0000000..5404b74 --- /dev/null +++ b/dev/check_fluids_can_be_called_from_perl.pl @@ -0,0 +1,145 @@ +#!/usr/bin/perl +use strict; +use warnings; +use Inline Python => <<'END_PYTHON'; +import fluids +END_PYTHON +use Time::HiRes qw(time); +use POSIX qw(ceil floor); + +# Function to test basic fluids functionality +sub test_fluids { + eval { + # Test basic module import and version + my $version = Inline::Python::py_eval('fluids.__version__', 0); + print "✓ Successfully imported fluids\n"; + print "✓ Fluids version: $version\n"; + + # Test basic Reynolds number calculation + # Create kwargs dictionary in Python + Inline::Python::py_eval('kwargs = {"V": 2.5, "D": 0.1, "rho": 1000, "mu": 0.001}'); + my $Re = Inline::Python::py_eval('fluids.Reynolds(**kwargs)', 0); + print "✓ Reynolds number calculation successful: $Re\n"; + die "Invalid Reynolds number" unless $Re > 0; + + # Test friction factor calculation + Inline::Python::py_eval('kwargs = {"Re": 1e5, "eD": 0.0001}'); + my $fd = Inline::Python::py_eval('fluids.friction_factor(**kwargs)', 0); + print "✓ Friction factor calculation successful: $fd\n"; + die "Invalid friction factor" unless ($fd > 0 && $fd < 1); + + print "\nAll basic tests completed successfully!\n"; + }; + if ($@) { + print "Error occurred: $@\n"; + die $@; + } +} + +# Function to test atmosphere calculations +sub test_atmosphere { + eval { + # Test ATMOSPHERE_1976 class + Inline::Python::py_eval('atm = fluids.ATMOSPHERE_1976(Z=5000)'); + + print "\nTesting atmosphere at 5000m elevation:\n"; + printf "✓ Temperature: %.4f\n", Inline::Python::py_eval('atm.T', 0); + printf "✓ Pressure: %.4f\n", Inline::Python::py_eval('atm.P', 0); + printf "✓ Density: %.6f\n", Inline::Python::py_eval('atm.rho', 0); + + # Test derived properties + printf "✓ Gravity: %.6f\n", Inline::Python::py_eval('atm.g', 0); + printf "✓ Viscosity: %.6e\n", Inline::Python::py_eval('atm.mu', 0); + printf "✓ Thermal conductivity: %.6f\n", Inline::Python::py_eval('atm.k', 0); + printf "✓ Sonic velocity: %.4f\n", Inline::Python::py_eval('atm.v_sonic', 0); + + # Test static methods + my $g_high = Inline::Python::py_eval('fluids.ATMOSPHERE_1976.gravity(Z=1E5)', 0); + printf "✓ High altitude gravity: %.6f\n", $g_high; + + my $v_sonic = Inline::Python::py_eval('fluids.ATMOSPHERE_1976.sonic_velocity(T=300)', 0); + printf "✓ Sonic velocity at 300K: %.4f\n", $v_sonic; + }; + if ($@) { + print "Error in atmosphere tests: $@\n"; + die $@; + } +} + +# Function to test tank calculations + + +# Function to test tank calculations +sub test_tank { + eval { + # Test basic tank creation + Inline::Python::py_eval(<<'END_PYTHON'); +T1 = fluids.TANK(V=10, L_over_D=0.7, sideB='conical', horizontal=False) +END_PYTHON + print "\nTesting tank calculations:\n"; + printf "✓ Tank length: %.6f\n", Inline::Python::py_eval('T1.L', 0); + printf "✓ Tank diameter: %.6f\n", Inline::Python::py_eval('T1.D', 0); + + # Test ellipsoidal tank + Inline::Python::py_eval(<<'END_PYTHON'); +tank_ellip = fluids.TANK(D=10, V=500, horizontal=False, + sideA='ellipsoidal', sideB='ellipsoidal', + sideA_a=1, sideB_a=1) +END_PYTHON + printf "✓ Ellipsoidal tank L: %.6f\n", Inline::Python::py_eval('tank_ellip.L', 0); + + # Test torispherical tank + Inline::Python::py_eval(<<'END_PYTHON'); +DIN = fluids.TANK(L=3, D=5, horizontal=False, + sideA='torispherical', sideB='torispherical', + sideA_f=1, sideA_k=0.1, sideB_f=1, sideB_k=0.1) +END_PYTHON + printf "✓ Tank max height: %.6f\n", Inline::Python::py_eval('DIN.h_max', 0); + printf "✓ Height at V=40: %.6f\n", Inline::Python::py_eval('DIN.h_from_V(40)', 0); + printf "✓ Volume at h=4.1: %.5f\n", Inline::Python::py_eval('DIN.V_from_h(4.1)', 0); + printf "✓ Surface area at h=2.1: %.5f\n", Inline::Python::py_eval('DIN.SA_from_h(2.1)', 0); + }; + if ($@) { + print "Error in tank tests: $@\n"; + die $@; + } +} +# Function to benchmark fluids operations + + +# Function to benchmark fluids operations +sub benchmark_fluids { + print "\nRunning benchmarks:\n"; + + # Benchmark friction factor calculation + print "\nBenchmarking friction_factor:\n"; + my $start_time = time(); + Inline::Python::py_eval(<<'END_PYTHON'); +for i in range(10000): + fluids.friction_factor(Re=1e5, eD=0.0001) +END_PYTHON + my $duration = time() - $start_time; + printf "Time for 10000 friction_factor calls: %.6f seconds\n", $duration; + printf "Average time per call: %.6f seconds\n", $duration/10000; + + # Benchmark tank creation + print "\nBenchmarking TANK creation:\n"; + $start_time = time(); + Inline::Python::py_eval(<<'END_PYTHON'); +for i in range(1000): + fluids.TANK(L=3, D=5, horizontal=False, + sideA='torispherical', sideB='torispherical', + sideA_f=1, sideA_k=0.1, sideB_f=1, sideB_k=0.1) +END_PYTHON + $duration = time() - $start_time; + printf "Time for 1000 TANK creations: %.6f seconds\n", $duration; + printf "Average time per creation: %.6f seconds\n", $duration/1000; +} + +# Run all tests +print "Running fluids tests from Perl...\n"; +test_fluids(); +test_atmosphere(); +test_tank(); +benchmark_fluids(); +print "\nAll tests completed!\n"; diff --git a/dev/check_fluids_can_be_called_from_ruby.rb b/dev/check_fluids_can_be_called_from_ruby.rb new file mode 100644 index 0000000..1aeac1e --- /dev/null +++ b/dev/check_fluids_can_be_called_from_ruby.rb @@ -0,0 +1,199 @@ +require 'pycall' + +# Import fluids +Fluids = PyCall.import_module('fluids') + +def test_fluids + begin + # 1. Test basic module import + puts "✓ Successfully imported fluids" + puts "✓ Fluids version: #{Fluids.__version__}" + + # 2. Test basic Reynolds number calculation + re = Fluids.Reynolds(V: 2.5, D: 0.1, rho: 1000, mu: 0.001) + puts "✓ Reynolds number calculation successful: #{re}" + raise unless re > 0 + + # 3. Test friction factor calculation + fd = Fluids.friction_factor(Re: 1e5, eD: 0.0001) + puts "✓ Friction factor calculation successful: #{fd}" + raise unless (0 < fd && fd < 1) + + puts "\nAll basic tests completed successfully!" + rescue => e + puts "Error occurred: #{e}" + raise + end +end + +def test_atmosphere + begin + # Test ATMOSPHERE_1976 class + atm = Fluids.ATMOSPHERE_1976.new(Z: 5000) + + puts "\nTesting atmosphere at 5000m elevation:" + puts "✓ Temperature: #{PyCall.getattr(atm, :T).round(4)}" + puts "✓ Pressure: #{PyCall.getattr(atm, :P).round(4)}" + puts "✓ Density: #{PyCall.getattr(atm, :rho).round(6)}" + + # Test derived properties + puts "✓ Gravity: #{PyCall.getattr(atm, :g).round(6)}" + puts format("✓ Viscosity: %.6e", PyCall.getattr(atm, :mu)) + puts "✓ Thermal conductivity: #{PyCall.getattr(atm, :k).round(6)}" + puts "✓ Sonic velocity: #{PyCall.getattr(atm, :v_sonic).round(4)}" + + # Test static methods + g_high = Fluids.ATMOSPHERE_1976.gravity(Z: 1E5) + puts "✓ High altitude gravity: #{g_high.round(6)}" + + v_sonic = Fluids.ATMOSPHERE_1976.sonic_velocity(T: 300) + puts "✓ Sonic velocity at 300K: #{v_sonic.round(4)}" + + mu_400 = Fluids.ATMOSPHERE_1976.viscosity(T: 400) + puts format("✓ Viscosity at 400K: %.6e", mu_400) + + k_400 = Fluids.ATMOSPHERE_1976.thermal_conductivity(T: 400) + puts "✓ Thermal conductivity at 400K: #{k_400.round(6)}" + rescue => e + puts "Error in atmosphere tests: #{e}" + raise + end +end + +def test_tank + begin + # Test basic tank creation + t1 = Fluids.TANK.new(V: 10, L_over_D: 0.7, sideB: "conical", horizontal: false) + puts "\nTesting tank calculations:" + puts "✓ Tank length: #{t1.L.round(6)}" + puts "✓ Tank diameter: #{t1.D.round(6)}" + + # Test ellipsoidal tank + tank_ellip = Fluids.TANK.new(D: 10, V: 500, horizontal: false, + sideA: "ellipsoidal", sideB: "ellipsoidal", + sideA_a: 1, sideB_a: 1) + puts "✓ Ellipsoidal tank L: #{tank_ellip.L.round(6)}" + + # Test torispherical tank + din = Fluids.TANK.new(L: 3, D: 5, horizontal: false, + sideA: "torispherical", sideB: "torispherical", + sideA_f: 1, sideA_k: 0.1, sideB_f: 1, sideB_k: 0.1) + + puts "✓ Tank representation: #{din}" + puts "✓ Tank max height: #{din.h_max.round(6)}" + puts "✓ Height at V=40: #{din.h_from_V(40).round(6)}" + puts "✓ Volume at h=4.1: #{din.V_from_h(4.1).round(5)}" + puts "✓ Surface area at h=2.1: #{din.SA_from_h(2.1).round(5)}" + rescue => e + puts "Error in tank tests: #{e}" + raise + end +end + +def test_reynolds + begin + puts "\nTesting Reynolds number calculations:" + + # Test with density and viscosity + re1 = Fluids.Reynolds(V: 2.5, D: 0.25, rho: 1.1613, mu: 1.9E-5) + puts "✓ Re (with rho, mu): #{re1.round(4)}" + raise unless (re1 - 38200.6579).abs < 0.1 + + # Test with kinematic viscosity + re2 = Fluids.Reynolds(V: 2.5, D: 0.25, nu: 1.636e-05) + puts "✓ Re (with nu): #{re2.round(4)}" + raise unless (re2 - 38202.934).abs < 0.1 + rescue => e + puts "Error in Reynolds tests: #{e}" + raise + end +end + +def test_psd + begin + puts "\nTesting particle size distribution functionality:" + + # Create a discrete PSD + ds = [240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532] + numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1] + + psd = Fluids.particle_size_distribution.ParticleSizeDistribution.new( + ds: ds, + fractions: numbers, + order: 0 + ) + puts "✓ Created discrete PSD" + + # Test mean sizes + d21 = psd.mean_size(2, 1) + puts "✓ Size-weighted mean diameter: #{d21.round(4)}" + raise unless (d21 - 1857.788).abs < 0.1 + + d10 = psd.mean_size(1, 0) + puts "✓ Arithmetic mean diameter: #{d10.round(4)}" + raise unless (d10 - 1459.372).abs < 0.1 + + # Test percentile calculations + d10_percentile = psd.dn(0.1) + d90_percentile = psd.dn(0.9) + puts "✓ D10: #{d10_percentile.round(4)}" + puts "✓ D90: #{d90_percentile.round(4)}" + + # Test probability functions + pdf_val = psd.pdf(1000) + cdf_val = psd.cdf(5000) + puts "✓ PDF at 1000: %.4e" % pdf_val + puts "✓ CDF at 5000: #{cdf_val.round(6)}" + + # Test lognormal distribution + psd_log = Fluids.particle_size_distribution.PSDLognormal.new(s: 0.5, d_characteristic: 5E-6) + puts "✓ Created lognormal PSD" + + vssa = psd_log.vssa + puts "✓ Volume specific surface area: #{vssa.round(2)}" + + span = psd_log.dn(0.9) - psd_log.dn(0.1) + puts "✓ Span: %.4e" % span + + ratio_7525 = psd_log.dn(0.75)/psd_log.dn(0.25) + puts "✓ D75/D25 ratio: #{ratio_7525.round(6)}" + rescue => e + puts "Error in PSD tests: #{e}" + raise + end +end + +def benchmark_fluids + puts "\nRunning benchmarks:" + + # Benchmark friction factor calculation + puts "\nBenchmarking friction_factor:" + t1 = Time.now + 1_000_000.times do + Fluids.friction_factor(Re: 1e5, eD: 0.0001) + end + elapsed = Time.now - t1 + puts "Time for 1e6 friction_factor calls: #{elapsed.round(6)} seconds" + puts "Average time per call: #{(elapsed/1_000_000).round(6)} seconds" + + # Benchmark tank creation + puts "\nBenchmarking TANK creation:" + t2 = Time.now + 1000.times do + Fluids.TANK.new(L: 3, D: 5, horizontal: false, + sideA: "torispherical", sideB: "torispherical", + sideA_f: 1, sideA_k: 0.1, sideB_f: 1, sideB_k: 0.1) + end + tank_elapsed = Time.now - t2 + puts "Average time per creation: #{(tank_elapsed/1000).round(6)} seconds" +end + +# Run all tests +puts "Running fluids tests from Ruby..." +test_fluids +test_atmosphere +test_tank +test_reynolds +test_psd +benchmark_fluids +puts "\nAll tests completed!" diff --git a/dev/cpp/CMakeLists.txt b/dev/cpp/CMakeLists.txt new file mode 100644 index 0000000..a20de96 --- /dev/null +++ b/dev/cpp/CMakeLists.txt @@ -0,0 +1,13 @@ +cmake_minimum_required(VERSION 3.15) +project(fluids_pybind11_test) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +find_package(Python COMPONENTS Interpreter Development REQUIRED) +find_package(pybind11 CONFIG REQUIRED) +message(STATUS "Python_EXECUTABLE: ${Python_EXECUTABLE}") +message(STATUS "Python_INCLUDE_DIRS: ${Python_INCLUDE_DIRS}") +message(STATUS "Python_LIBRARIES: ${Python_LIBRARIES}") +add_executable(fluids_test main.cpp) +target_link_libraries(fluids_test PRIVATE pybind11::embed) diff --git a/dev/cpp/main.cpp b/dev/cpp/main.cpp new file mode 100644 index 0000000..77ee1b6 --- /dev/null +++ b/dev/cpp/main.cpp @@ -0,0 +1,286 @@ +#include +#include // Added this for interpreter management +#include +#include +#include +#include +#include + +namespace py = pybind11; +PYBIND11_EMBEDDED_MODULE(dummy, m) {} // Needed for embedding + +class FluidsTester { +private: + py::scoped_interpreter guard; // RAII management of the Python interpreter + py::module_ fluids; + +public: + FluidsTester() { + try { + // The interpreter is automatically initialized by scoped_interpreter + fluids = py::module_::import("fluids"); + std::cout << "✓ Successfully imported fluids\n"; + std::cout << "✓ Fluids version: " << fluids.attr("__version__").cast() << "\n"; + } catch (const std::exception& e) { + std::cerr << "Error initializing Python: " << e.what() << std::endl; + throw; + } + } + + void test_fluids() { + try { + // Test basic Reynolds number calculation + auto Re = fluids.attr("Reynolds")( + py::arg("V") = 2.5, + py::arg("D") = 0.1, + py::arg("rho") = 1000, + py::arg("mu") = 0.001 + ).cast(); + std::cout << "✓ Reynolds number calculation successful: " << Re << "\n"; + assert(Re > 0); + + // Test friction factor calculation + auto fd = fluids.attr("friction_factor")( + py::arg("Re") = 1e5, + py::arg("eD") = 0.0001 + ).cast(); + std::cout << "✓ Friction factor calculation successful: " << fd << "\n"; + assert(fd > 0 && fd < 1); + + std::cout << "\nAll basic tests completed successfully!\n"; + } catch (const std::exception& e) { + std::cerr << "Error in fluids tests: " << e.what() << std::endl; + throw; + } + } + + void test_atmosphere() { + try { + // Test ATMOSPHERE_1976 class + auto atm = fluids.attr("ATMOSPHERE_1976")(py::arg("Z") = 5000); + + std::cout << "\nTesting atmosphere at 5000m elevation:\n"; + std::cout << "✓ Temperature: " << std::fixed << std::setprecision(4) + << atm.attr("T").cast() << "\n"; + std::cout << "✓ Pressure: " << atm.attr("P").cast() << "\n"; + std::cout << "✓ Density: " << std::setprecision(6) + << atm.attr("rho").cast() << "\n"; + + // Test derived properties + std::cout << "✓ Gravity: " << atm.attr("g").cast() << "\n"; + std::cout << "✓ Viscosity: " << std::scientific + << atm.attr("mu").cast() << "\n"; + std::cout << "✓ Thermal conductivity: " << std::fixed + << atm.attr("k").cast() << "\n"; + std::cout << "✓ Sonic velocity: " << std::setprecision(4) + << atm.attr("v_sonic").cast() << "\n"; + + // Test static methods + auto atm_class = fluids.attr("ATMOSPHERE_1976"); + auto g_high = atm_class.attr("gravity")(py::arg("Z") = 1E5).cast(); + std::cout << "✓ High altitude gravity: " << std::setprecision(6) + << g_high << "\n"; + + } catch (const std::exception& e) { + std::cerr << "Error in atmosphere tests: " << e.what() << std::endl; + throw; + } + } + + void test_tank() { + try { + // Test basic tank creation + auto T1 = fluids.attr("TANK")( + py::arg("V") = 10, + py::arg("L_over_D") = 0.7, + py::arg("sideB") = "conical", + py::arg("horizontal") = false + ); + + std::cout << "\nTesting tank calculations:\n"; + std::cout << "✓ Tank length: " << std::fixed << std::setprecision(6) + << T1.attr("L").cast() << "\n"; + std::cout << "✓ Tank diameter: " << T1.attr("D").cast() << "\n"; + + // Test torispherical tank + auto DIN = fluids.attr("TANK")( + py::arg("L") = 3, + py::arg("D") = 5, + py::arg("horizontal") = false, + py::arg("sideA") = "torispherical", + py::arg("sideB") = "torispherical", + py::arg("sideA_f") = 1, + py::arg("sideA_k") = 0.1, + py::arg("sideB_f") = 1, + py::arg("sideB_k") = 0.1 + ); + + std::cout << "✓ Tank max height: " << DIN.attr("h_max").cast() << "\n"; + std::cout << "✓ Height at V=40: " << DIN.attr("h_from_V")(40).cast() << "\n"; + std::cout << "✓ Volume at h=4.1: " << std::setprecision(5) + << DIN.attr("V_from_h")(4.1).cast() << "\n"; + + } catch (const std::exception& e) { + std::cerr << "Error in tank tests: " << e.what() << std::endl; + throw; + } + } +void test_reynolds() { + try { + std::cout << "\nTesting Reynolds number calculations:\n"; + + // Test with density and viscosity + auto Re1 = fluids.attr("Reynolds")( + py::arg("V") = 2.5, + py::arg("D") = 0.25, + py::arg("rho") = 1.1613, + py::arg("mu") = 1.9E-5 + ).cast(); + std::cout << "✓ Re (with rho, mu): " << std::fixed << std::setprecision(4) << Re1 << "\n"; + assert(std::abs(Re1 - 38200.6579) < 0.1); + + // Test with kinematic viscosity + auto Re2 = fluids.attr("Reynolds")( + py::arg("V") = 2.5, + py::arg("D") = 0.25, + py::arg("nu") = 1.636e-05 + ).cast(); + std::cout << "✓ Re (with nu): " << Re2 << "\n"; + assert(std::abs(Re2 - 38202.934) < 0.1); + + } catch (const std::exception& e) { + std::cerr << "Error in Reynolds tests: " << e.what() << std::endl; + throw; + } +} + +void test_psd() { + try { + std::cout << "\nTesting particle size distribution functionality:\n"; + + // Create vectors for discrete PSD + std::vector ds = {240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, + 2141, 2676, 3345, 4181, 5226, 6532}; + std::vector numbers = {65, 119, 232, 410, 629, 849, 990, 981, 825, + 579, 297, 111, 21, 1}; + + // Create a discrete PSD + auto psd = fluids.attr("particle_size_distribution").attr("ParticleSizeDistribution")( + py::arg("ds") = ds, + py::arg("fractions") = numbers, + py::arg("order") = 0 + ); + std::cout << "✓ Created discrete PSD\n"; + + // Test mean sizes + auto d21 = psd.attr("mean_size")(2, 1).cast(); + std::cout << "✓ Size-weighted mean diameter: " << std::fixed + << std::setprecision(4) << d21 << "\n"; + assert(std::abs(d21 - 1857.788) < 0.1); + + auto d10 = psd.attr("mean_size")(1, 0).cast(); + std::cout << "✓ Arithmetic mean diameter: " << d10 << "\n"; + assert(std::abs(d10 - 1459.372) < 0.1); + + // Test percentile calculations + auto d10_percentile = psd.attr("dn")(0.1).cast(); + auto d90_percentile = psd.attr("dn")(0.9).cast(); + std::cout << "✓ D10: " << d10_percentile << "\n"; + std::cout << "✓ D90: " << d90_percentile << "\n"; + + // Test probability functions + auto pdf_val = psd.attr("pdf")(1000).cast(); + auto cdf_val = psd.attr("cdf")(5000).cast(); + std::cout << "✓ PDF at 1000: " << std::scientific << pdf_val << "\n"; + std::cout << "✓ CDF at 5000: " << std::fixed << std::setprecision(6) << cdf_val << "\n"; + + // Test lognormal distribution + auto psd_log = fluids.attr("particle_size_distribution").attr("PSDLognormal")( + py::arg("s") = 0.5, + py::arg("d_characteristic") = 5E-6 + ); + std::cout << "✓ Created lognormal PSD\n"; + + auto vssa = psd_log.attr("vssa").cast(); + std::cout << "✓ Volume specific surface area: " << std::setprecision(2) + << vssa << "\n"; + + auto span = psd_log.attr("dn")(0.9).cast() - + psd_log.attr("dn")(0.1).cast(); + std::cout << "✓ Span: " << std::scientific << span << "\n"; + + auto ratio_7525 = (psd_log.attr("dn")(0.75).cast() / + psd_log.attr("dn")(0.25).cast()); + std::cout << "✓ D75/D25 ratio: " << std::fixed << std::setprecision(6) + << ratio_7525 << "\n"; + + } catch (const std::exception& e) { + std::cerr << "Error in PSD tests: " << e.what() << std::endl; + throw; + } +} + void benchmark_fluids() { + std::cout << "\nRunning benchmarks:\n"; + + // Benchmark friction factor calculation + std::cout << "\nBenchmarking friction_factor:\n"; + auto start = std::chrono::high_resolution_clock::now(); + + for(int i = 0; i < 1000000; i++) { + fluids.attr("friction_factor")( + py::arg("Re") = 1e5, + py::arg("eD") = 0.0001 + ); + } + + auto end = std::chrono::high_resolution_clock::now(); + auto duration = std::chrono::duration_cast(end - start); + + std::cout << "Time for 1e6 friction_factor calls: " + << duration.count() << " microseconds\n"; + std::cout << "Average time per call: " + << duration.count() / 1000000.0 << " microseconds\n"; + + // Benchmark tank creation + std::cout << "\nBenchmarking TANK creation:\n"; + start = std::chrono::high_resolution_clock::now(); + + for(int i = 0; i < 1000; i++) { + fluids.attr("TANK")( + py::arg("L") = 3, + py::arg("D") = 5, + py::arg("horizontal") = false, + py::arg("sideA") = "torispherical", + py::arg("sideB") = "torispherical", + py::arg("sideA_f") = 1, + py::arg("sideA_k") = 0.1, + py::arg("sideB_f") = 1, + py::arg("sideB_k") = 0.1 + ); + } + + end = std::chrono::high_resolution_clock::now(); + duration = std::chrono::duration_cast(end - start); + + std::cout << "Average time per tank creation: " + << duration.count() / 1000.0 << " microseconds\n"; + } +}; + +int main() { + try { + std::cout << "Running fluids tests from C++...\n"; + FluidsTester tester; + tester.test_fluids(); + tester.test_atmosphere(); + tester.test_tank(); + tester.test_reynolds(); + tester.test_psd(); + tester.benchmark_fluids(); + std::cout << "\nAll tests completed!\n"; + return 0; + } catch (const std::exception& e) { + std::cerr << "Fatal error: " << e.what() << std::endl; + return 1; + } +} diff --git a/dev/haskell-fluids-test/Main.hs b/dev/haskell-fluids-test/Main.hs deleted file mode 100644 index 5d3e94d..0000000 --- a/dev/haskell-fluids-test/Main.hs +++ /dev/null @@ -1,26 +0,0 @@ --- module Main where --- --- main :: IO () --- main = putStrLn "Hello, Haskell!" -{-# LANGUAGE OverloadedStrings #-} - -module Main where - -import qualified CPython as Py -import qualified CPython.Simple as PySim - -main :: IO () -main = do - putStrLn "Testing fluids from Haskell..." - -- Initialize Python - Py.initialize - - -- Basic test - re <- PySim.call "fluids" "Reynolds" - [PySim.arg (2.5 :: Double), PySim.arg (0.1 :: Double), - PySim.arg (1000.0 :: Double), PySim.arg (0.001 :: Double)] [] - - putStrLn $ "Reynolds number calculation: " ++ show (re :: Double) - - -- Clean up - Py.finalize diff --git a/dev/haskell-fluids-test/app/Main.hs b/dev/haskell-fluids-test/app/Main.hs deleted file mode 100644 index 52b13c3..0000000 --- a/dev/haskell-fluids-test/app/Main.hs +++ /dev/null @@ -1,75 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE ScopedTypeVariables #-} - -module Main where - -import qualified CPython as Py -import qualified CPython.Simple as PySim -import Data.Time.Clock.System (getSystemTime, systemSeconds) -import Control.Exception (handle) -import qualified CPython.Types.Exception as PyExc -import Text.Printf (printf) -import Data.Text (Text) - --- Helper function for timing -time :: IO a -> IO (Double, a) -time action = do - start <- getSystemTime - result <- action - end <- getSystemTime - let diff = fromIntegral (systemSeconds end - systemSeconds start) - return (diff, result) - --- Basic fluids test -testFluids :: IO () -testFluids = do - putStrLn "Testing basic fluids functionality..." - - -- Test version - version <- (PySim.getAttribute "fluids" "__version__" :: IO Text) - putStrLn $ "✓ Fluids version: " ++ show version - - -- Test Reynolds number - re <- (PySim.call "fluids" "Reynolds" - [PySim.arg (2.5 :: Double), PySim.arg (0.1 :: Double), - PySim.arg (1000.0 :: Double), PySim.arg (0.001 :: Double)] [] :: IO Double) - putStrLn $ "✓ Reynolds number calculation: " ++ show re - - -- Test friction factor - fd <- (PySim.call "fluids" "friction_factor" - [PySim.arg (1e5 :: Double), PySim.arg (0.0001 :: Double)] [] :: IO Double) - putStrLn $ "✓ Friction factor calculation: " ++ show fd - - putStrLn "Basic tests completed successfully!\n" - - --- Benchmark fluids functions -benchmarkFluids :: IO () -benchmarkFluids = do - putStrLn "Running benchmarks..." - - -- Benchmark friction factor - putStrLn "\nBenchmarking friction_factor:" - (t1, _) <- time $ do - sequence $ replicate 1000 $ - (PySim.call "fluids" "friction_factor" - [PySim.arg (1e5 :: Double), PySim.arg (0.0001 :: Double)] [] :: IO Double) - putStrLn $ printf "Time for 1000 friction_factor calls: %.6f seconds" t1 - putStrLn $ printf "Average time per call: %.6f seconds" (t1 / 1000) - --- Main function to run all tests -main :: IO () -main = handle pyExceptionHandler $ do - putStrLn "Running fluids tests from Haskell..." - Py.initialize - - testFluids - benchmarkFluids - - Py.finalize - putStrLn "All tests completed!" - where - pyExceptionHandler :: PyExc.Exception -> IO () - pyExceptionHandler e = do - putStrLn $ "Python error occurred: " ++ show e - Py.finalize diff --git a/dev/haskell-fluids-test/CHANGELOG.md b/dev/haskell/CHANGELOG.md similarity index 100% rename from dev/haskell-fluids-test/CHANGELOG.md rename to dev/haskell/CHANGELOG.md diff --git a/dev/haskell/app/Main.hs b/dev/haskell/app/Main.hs new file mode 100644 index 0000000..2dce76e --- /dev/null +++ b/dev/haskell/app/Main.hs @@ -0,0 +1,165 @@ +{-# LANGUAGE OverloadedStrings #-} +module Main where + +import qualified CPython as Py +import qualified CPython.Types.Module as Py +import qualified CPython.Protocols.Object as Py +import qualified CPython.Protocols.Number as PyNum +import qualified CPython.Types.Dictionary as PyDict +import qualified CPython.Types.Tuple as PyTuple +import qualified CPython.Types.Unicode as PyUnicode +import qualified CPython.Types.Float as PyFloat +import qualified CPython.Types.Exception as PyExc +import Control.Exception (handle) +import Text.Printf (printf) +import Data.Maybe (fromMaybe) +import qualified Data.Text as T +import System.Clock + +testBasics :: IO () +testBasics = do + putStrLn "Testing basic fluids functionality..." + + -- Import fluids module + fluidsModule <- Py.importModule (T.pack "fluids") + + -- Get version + versionName <- PyUnicode.toUnicode (T.pack "__version__") + versionObj <- Py.getAttribute fluidsModule versionName + Just versionStr <- Py.cast versionObj + version <- PyUnicode.fromUnicode versionStr + putStrLn $ "✓ Fluids version: " ++ T.unpack version + + -- Test Reynolds number calculation + let calcReynolds = do + -- Get Reynolds function + reynoldsName <- PyUnicode.toUnicode (T.pack "Reynolds") + reynolds <- Py.getAttribute fluidsModule reynoldsName + + -- Create arguments + v <- PyFloat.toFloat 2.5 + d <- PyFloat.toFloat 0.1 + rho <- PyFloat.toFloat 1000.0 + mu <- PyFloat.toFloat 0.001 + args <- PyTuple.toTuple [Py.toObject v, Py.toObject d, Py.toObject rho, Py.toObject mu] + kwargs <- PyDict.new + + -- Call function + result <- Py.call reynolds args kwargs + x <- PyNum.castToNumber result + let num = fromMaybe (error "Could not convert Reynolds number") x + PyFloat.fromFloat =<< PyNum.toFloat num + + re <- calcReynolds + putStrLn $ printf "✓ Reynolds number calculation: %.1f" re + + -- Test friction factor + let calcFriction = do + -- Get friction_factor function + frictionName <- PyUnicode.toUnicode (T.pack "friction_factor") + friction <- Py.getAttribute fluidsModule frictionName + + -- Create arguments + re <- PyFloat.toFloat 1e5 + ed <- PyFloat.toFloat 0.0001 + args <- PyTuple.toTuple [Py.toObject re, Py.toObject ed] + kwargs <- PyDict.new + + -- Call function + result <- Py.call friction args kwargs + x <- PyNum.castToNumber result + let num = fromMaybe (error "Could not convert friction factor") x + PyFloat.fromFloat =<< PyNum.toFloat num + + fd <- calcFriction + putStrLn $ printf "✓ Friction factor calculation: %.6f" fd + putStrLn "Basic tests completed successfully!\n" + +testAtmosphere :: IO () +testAtmosphere = do + putStrLn "\nTesting atmosphere at 5000m elevation:" + + -- Import fluids module + fluidsModule <- Py.importModule (T.pack "fluids") + + -- Create argument for constructor + zArg <- PyFloat.toFloat 5000.0 + args <- PyTuple.toTuple [Py.toObject zArg] + kwargs <- PyDict.new + + -- Get ATMOSPHERE_1976 class and create instance + atmosClass <- PyUnicode.toUnicode (T.pack "ATMOSPHERE_1976") >>= Py.getAttribute fluidsModule + atm <- Py.call atmosClass args kwargs + + -- Get and print properties + let getProperty name = do + nameObj <- PyUnicode.toUnicode name + prop <- Py.getAttribute atm nameObj + x <- PyNum.castToNumber prop + let num = fromMaybe (error $ "Could not get " ++ T.unpack name ++ " as number") x + PyFloat.fromFloat =<< PyNum.toFloat num + + temp <- getProperty (T.pack "T") + pressure <- getProperty (T.pack "P") + density <- getProperty (T.pack "rho") + gravity <- getProperty (T.pack "g") + viscosity <- getProperty (T.pack "mu") + conductivity <- getProperty (T.pack "k") + sonicVel <- getProperty (T.pack "v_sonic") + + putStrLn $ printf "✓ Temperature: %.4f" temp + putStrLn $ printf "✓ Pressure: %.4f" pressure + putStrLn $ printf "✓ Density: %.6f" density + putStrLn $ printf "✓ Gravity: %.6f" gravity + putStrLn $ printf "✓ Viscosity: %.6e" viscosity + putStrLn $ printf "✓ Thermal conductivity: %.6f" conductivity + putStrLn $ printf "✓ Sonic velocity: %.4f" sonicVel + +benchmarkFluids :: IO () +benchmarkFluids = do + putStrLn "\nRunning benchmarks:" + + -- Import fluids module + fluidsModule <- Py.importModule (T.pack "fluids") + + -- Get friction_factor function + frictionName <- PyUnicode.toUnicode (T.pack "friction_factor") + friction <- Py.getAttribute fluidsModule frictionName + + -- Prepare arguments that will be reused + re <- PyFloat.toFloat 1e5 + ed <- PyFloat.toFloat 0.0001 + args <- PyTuple.toTuple [Py.toObject re, Py.toObject ed] + kwargs <- PyDict.new + + -- Time the operations + putStrLn "\nBenchmarking friction_factor:" + start <- getTime Monotonic + + -- Do 1,000,000 iterations like in Julia version + sequence_ $ replicate 1000000 $ do + result <- Py.call friction args kwargs + x <- PyNum.castToNumber result + let num = fromMaybe (error "Could not convert friction factor") x + PyFloat.fromFloat =<< PyNum.toFloat num + + end <- getTime Monotonic + let diff = (fromIntegral (toNanoSecs (diffTimeSpec end start)) :: Double) / 1000000000 + + putStrLn $ printf "Time for 1e6 friction_factor calls: %.6f seconds" diff + putStrLn $ printf "Average time per call: %.6f seconds" (diff / 1000000) + +main :: IO () +main = handle pyExceptionHandler $ do + putStrLn "Running fluids tests from Haskell..." + Py.initialize + testBasics + testAtmosphere + benchmarkFluids + Py.finalize + putStrLn "\nAll tests completed!" + where + pyExceptionHandler :: PyExc.Exception -> IO () + pyExceptionHandler e = do + putStrLn $ "Python error occurred: " ++ show e + Py.finalize diff --git a/dev/haskell-fluids-test/haskell-fluids-test.cabal b/dev/haskell/haskell-fluids-test.cabal similarity index 92% rename from dev/haskell-fluids-test/haskell-fluids-test.cabal rename to dev/haskell/haskell-fluids-test.cabal index 284f07d..6a210f2 100644 --- a/dev/haskell-fluids-test/haskell-fluids-test.cabal +++ b/dev/haskell/haskell-fluids-test.cabal @@ -26,6 +26,7 @@ executable haskell-fluids-test build-depends: base ^>=4.15.1.0, cpython, time, - text + text, + clock hs-source-dirs: app default-language: Haskell2010 diff --git a/dev/nim/check_fluids_can_be_called_from_nim.nim b/dev/nim/check_fluids_can_be_called_from_nim.nim new file mode 100644 index 0000000..0058fec --- /dev/null +++ b/dev/nim/check_fluids_can_be_called_from_nim.nim @@ -0,0 +1,211 @@ +import nimpy +import times +import strformat +import math + +# Initialize Python and import fluids +let fluids = pyImport("fluids") + +proc test_fluids() = + try: + # Test basic module import + echo "Running fluids tests from Nim..." + echo "✓ Successfully imported fluids" + echo &"✓ Fluids version: {fluids.getAttr(\"__version__\").to(string)}" + + # Test basic Reynolds number calculation + let Re = fluids.Reynolds(V=2.5, D=0.1, rho=1000, mu=0.001).to(float) + echo &"✓ Reynolds number calculation successful: {Re}" + assert Re > 0 + + # Test friction factor calculation + let fd = fluids.friction_factor(Re=1e5, eD=0.0001).to(float) + echo &"✓ Friction factor calculation successful: {fd}" + assert 0 < fd and fd < 1 + + echo "\nAll basic tests completed successfully!" + + except: + echo "Error in fluids tests: ", getCurrentExceptionMsg() + raise + +proc test_atmosphere() = + try: + # Test ATMOSPHERE_1976 class + let atm = fluids.ATMOSPHERE_1976(Z=5000) + + echo "\nTesting atmosphere at 5000m elevation:" + echo &"✓ Temperature: {atm.T.to(float):.4f}" + echo &"✓ Pressure: {atm.P.to(float):.4f}" + echo &"✓ Density: {atm.rho.to(float):.6f}" + + # Test derived properties + echo &"✓ Gravity: {atm.g.to(float):.6f}" + echo &"✓ Viscosity: {atm.mu.to(float):.6e}" + echo &"✓ Thermal conductivity: {atm.k.to(float):.6f}" + echo &"✓ Sonic velocity: {atm.v_sonic.to(float):.4f}" + + # Test static methods + let g_high = fluids.ATMOSPHERE_1976.gravity(Z=1E5).to(float) + echo &"✓ High altitude gravity: {g_high:.6f}" + + let v_sonic = fluids.ATMOSPHERE_1976.sonic_velocity(T=300).to(float) + echo &"✓ Sonic velocity at 300K: {v_sonic:.4f}" + + let mu_400 = fluids.ATMOSPHERE_1976.viscosity(T=400).to(float) + echo &"✓ Viscosity at 400K: {mu_400:.6e}" + + let k_400 = fluids.ATMOSPHERE_1976.thermal_conductivity(T=400).to(float) + echo &"✓ Thermal conductivity at 400K: {k_400:.6f}" + + except: + echo "Error in atmosphere tests: ", getCurrentExceptionMsg() + raise + +proc test_tank() = + try: + # Test basic tank creation + let T1 = fluids.TANK(V=10, L_over_D=0.7, sideB="conical", horizontal=false) + echo "\nTesting tank calculations:" + echo &"✓ Tank length: {T1.L.to(float):.6f}" + echo &"✓ Tank diameter: {T1.D.to(float):.6f}" + + # Test ellipsoidal tank + let tank_ellip = fluids.TANK( + D=10, V=500, horizontal=false, + sideA="ellipsoidal", sideB="ellipsoidal", + sideA_a=1, sideB_a=1 + ) + echo &"✓ Ellipsoidal tank L: {tank_ellip.L.to(float):.6f}" + + # Test torispherical tank + let DIN = fluids.TANK( + L=3, D=5, horizontal=false, + sideA="torispherical", sideB="torispherical", + sideA_f=1, sideA_k=0.1, sideB_f=1, sideB_k=0.1 + ) + + echo &"✓ Tank representation: {$DIN}" + echo &"✓ Tank max height: {DIN.h_max.to(float):.6f}" + echo &"✓ Height at V=40: {DIN.h_from_V(40).to(float):.6f}" + echo &"✓ Volume at h=4.1: {DIN.V_from_h(4.1).to(float):.5f}" + echo &"✓ Surface area at h=2.1: {DIN.SA_from_h(2.1).to(float):.5f}" + + except: + echo "Error in tank tests: ", getCurrentExceptionMsg() + raise + +proc test_reynolds() = + try: + echo "\nTesting Reynolds number calculations:" + + # Test with density and viscosity + let Re1 = fluids.Reynolds(V=2.5, D=0.25, rho=1.1613, mu=1.9E-5).to(float) + echo &"✓ Re (with rho, mu): {Re1:.4f}" + assert abs(Re1 - 38200.6579) < 0.1 + + # Test with kinematic viscosity + let Re2 = fluids.Reynolds(V=2.5, D=0.25, nu=1.636e-05).to(float) + echo &"✓ Re (with nu): {Re2:.4f}" + assert abs(Re2 - 38202.934) < 0.1 + + except: + echo "Error in Reynolds tests: ", getCurrentExceptionMsg() + raise + +proc test_psd() = + try: + echo "\nTesting particle size distribution functionality:" + + # Create a discrete PSD + let ds = @[240.0, 360.0, 450.0, 562.5, 703.0, 878.0, 1097.0, 1371.0, + 1713.0, 2141.0, 2676.0, 3345.0, 4181.0, 5226.0, 6532.0] + let numbers = @[65.0, 119.0, 232.0, 410.0, 629.0, 849.0, 990.0, 981.0, + 825.0, 579.0, 297.0, 111.0, 21.0, 1.0] + + let psd = fluids.particle_size_distribution.ParticleSizeDistribution( + ds=ds, + fractions=numbers, + order=0 + ) + echo "✓ Created discrete PSD" + + # Test mean sizes + let d21 = psd.mean_size(2, 1).to(float) + echo &"✓ Size-weighted mean diameter: {d21:.4f}" + assert abs(d21 - 1857.788) < 0.1 + + let d10 = psd.mean_size(1, 0).to(float) + echo &"✓ Arithmetic mean diameter: {d10:.4f}" + assert abs(d10 - 1459.372) < 0.1 + + # Test percentile calculations + let d10_percentile = psd.dn(0.1).to(float) + let d90_percentile = psd.dn(0.9).to(float) + echo &"✓ D10: {d10_percentile:.4f}" + echo &"✓ D90: {d90_percentile:.4f}" + + # Test probability functions + let pdf_val = psd.pdf(1000).to(float) + let cdf_val = psd.cdf(5000).to(float) + echo &"✓ PDF at 1000: {pdf_val:.4e}" + echo &"✓ CDF at 5000: {cdf_val:.6f}" + + # Test lognormal distribution + let psd_log = fluids.particle_size_distribution.PSDLognormal(s=0.5, d_characteristic=5E-6) + echo "✓ Created lognormal PSD" + + let vssa = psd_log.vssa.to(float) + echo &"✓ Volume specific surface area: {vssa:.2f}" + + let span = psd_log.dn(0.9).to(float) - psd_log.dn(0.1).to(float) + echo &"✓ Span: {span:.4e}" + + let ratio_7525 = psd_log.dn(0.75).to(float) / psd_log.dn(0.25).to(float) + echo &"✓ D75/D25 ratio: {ratio_7525:.6f}" + + except: + echo "Error in PSD tests: ", getCurrentExceptionMsg() + raise + +proc benchmark_fluids() = + echo "\nRunning benchmarks:" + + # Benchmark friction factor calculation + echo "\nBenchmarking friction_factor:" + let start = epochTime() + + for i in 1..1_000_000: + discard fluids.friction_factor(Re=1e5, eD=0.0001) + + let duration = (epochTime() - start) * 1_000_000 # Convert to microseconds + echo &"Time for 1e6 friction_factor calls: {duration.int} microseconds" + echo &"Average time per call: {duration / 1_000_000:.6f} microseconds" + + # Benchmark tank creation + echo "\nBenchmarking TANK creation:" + let tank_start = epochTime() + + for i in 1..1000: + discard fluids.TANK( + L=3, D=5, horizontal=false, + sideA="torispherical", sideB="torispherical", + sideA_f=1, sideA_k=0.1, + sideB_f=1, sideB_k=0.1 + ) + + let tank_duration = (epochTime() - tank_start) * 1_000_000 # Convert to microseconds + echo &"Average time per tank creation: {tank_duration / 1000:.6f} microseconds" + +when isMainModule: + try: + test_fluids() + test_atmosphere() + test_tank() + test_reynolds() + test_psd() + benchmark_fluids() + echo "\nAll tests completed!" + except: + echo "Fatal error: ", getCurrentExceptionMsg() + quit(1)