diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/.DS_Store differ diff --git a/.vscode/.DS_Store b/.vscode/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/.vscode/.DS_Store differ diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json new file mode 100644 index 0000000..e1b4277 --- /dev/null +++ b/.vscode/c_cpp_properties.json @@ -0,0 +1,22 @@ +{ + "configurations": [ + { + "name": "Mac", + "includePath": [ + "${workspaceFolder}/**", + "/opt/homebrew/opt/python@3.12/Frameworks/Python.framework/Versions/3.12/include/python3.12" + ], + "defines": [], + "macFrameworkPath": [ + "/System/Library/Frameworks", + "/Library/Frameworks" + ], + "compilerPath": "/usr/bin/clang", + "cStandard": "c11", + "cppStandard": "c++17", + "intelliSenseMode": "macos-clang-x64" + } + ], + "version": 4 + } + \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..4c2091a --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,48 @@ +{ + // Path to the Python interpreter + "python.defaultInterpreterPath": "/opt/homebrew/opt/python@3.12/Frameworks/Python.framework/Versions/3.12/bin/python3", + + // Use Pylance as the language server + "python.languageServer": "Pylance", + + // Basic type checking mode + "python.analysis.typeCheckingMode": "basic", + + // Automatically search for import paths + "python.analysis.autoSearchPaths": true, + + // Use library code for types + "python.analysis.useLibraryCodeForTypes": true, + + // Code Runner extension settings + "code-runner.executorMap": { + "python": "python -u" + }, + "code-runner.runInTerminal": true, + "code-runner.saveFileBeforeRun": true, + + // Enable linting + "python.linting.enabled": true, + + // Use Ruff as the linter + "python.linting.ruffEnabled": true, + + // Disable Pylint + "python.linting.pylintEnabled": false, + + // Additional arguments for Ruff (empty in this case) + "python.linting.ruffArgs": [], + + // File associations to specify how certain files are treated + "files.associations": { + "libfoo.C": "cpp", + "utils.C": "cpp", + "helloPy.C": "cpp", + "ftcs.C": "cpp", + "heat.C": "cpp", + "dufrank.C": "cpp", + "args.C": "cpp", + "crankn.C": "cpp", + "exact.C": "cpp" + } +} diff --git a/Number.H b/Number.H deleted file mode 100644 index 0a1ddd9..0000000 --- a/Number.H +++ /dev/null @@ -1,121 +0,0 @@ -#include - -#include "Half.H" - -#ifndef FPTYPE -#define FPTYPE 2 -#endif - -// 0=half, 1=float, 2=double, 3=quad -#if FPTYPE == 0 -typedef half_float::half fpnumber; -#elif FPTYPE == 1 -typedef float fpnumber; -#elif FPTYPE == 2 -typedef double fpnumber; -#elif FPTYPE == 3 -typedef long double fpnumber; -#else -#error UNRECOGNIZED FPTYPE -#endif - -// Encapsulate a floating point number to easily track various metrics. -// Note: Doing this well is complicated by the mix of types in numerical -// statements (e.g. int*double) made explicit by modern compilers. -class Number { - public: - static void *operator new(std::size_t sz) { Number::nbytes += sz; return ::operator new(sz); }; - static void *operator new[](std::size_t sz) { Number::nbytes += sz; return ::operator new(sz); }; - static char const *counts_string() { static char counts[256]; - snprintf(counts, sizeof(counts), "Adds:%d, Mults:%d, Divs:%d, Bytes:%d", Number::nadds, Number::nmults, Number::ndivs, (int) Number::nbytes); - return counts; }; - inline Number() : x(0) {}; - inline Number(fpnumber _x) : x(_x) {}; -#if FPTYPE!=2 - inline Number(double _x) : x((fpnumber) _x) {}; -#endif - inline Number(int _x) : x((fpnumber) _x) {}; - - inline Number &operator=(const Number& rhs) { this->x = rhs.x; return *this; }; - inline operator fpnumber() const { return x; }; - inline operator int() const { return static_cast(x); }; -#if FPTYPE!=2 - inline operator double() const { return static_cast(x); }; -#endif - private: - static int nadds; - static int nmults; - static int ndivs; - static std::size_t nbytes; - fpnumber x; - - #pragma omp threadprivate(nadds) - #pragma omp threadprivate(nmults) - #pragma omp threadprivate(ndivs) - #pragma omp threadprivate(nbytes) - - // Various operators on Number w/mix of int - friend Number operator+(const Number& lhs, const Number& rhs) { Number::nadds++; return lhs.x + rhs.x; }; - friend Number operator+(const int& lhs, const Number& rhs) { Number::nadds++; return lhs + rhs.x; }; - friend Number operator+(const Number& lhs, const int& rhs) { Number::nadds++; return lhs.x + rhs; }; - friend Number operator+(const double& lhs, const Number& rhs) { Number::nadds++; return lhs + rhs.x; }; - friend Number operator+(const Number& lhs, const double& rhs) { Number::nadds++; return lhs.x + rhs; }; - friend Number operator+=(Number& lhs, const Number& rhs) { Number::nadds++; return lhs.x += rhs.x; }; - friend Number operator+=(Number& lhs, const double& rhs) { Number::nadds++; return lhs.x += rhs; }; - friend Number operator+=(double& lhs, const Number& rhs) { Number::nadds++; return lhs += rhs.x; }; - friend Number operator-(const Number& lhs, const Number& rhs) { Number::nadds++; return lhs.x - rhs.x; }; - friend Number operator-(const int& lhs, const Number& rhs) { Number::nadds++; return lhs - rhs.x; }; - friend Number operator-(const Number& lhs, const int& rhs) { Number::nadds++; return lhs.x - rhs; }; - friend Number operator-(const double& lhs, const Number& rhs) { Number::nadds++; return lhs - rhs.x; }; - friend Number operator-(const Number& lhs, const double& rhs) { Number::nadds++; return lhs.x - rhs; }; - friend Number operator-(const Number& rhs) { Number::nadds++; return -rhs.x; }; - friend Number operator-=(Number& lhs, const Number& rhs) { Number::nadds++; return lhs.x -= rhs.x; }; - friend Number operator*(const Number& lhs, const Number& rhs) { Number::nmults++; return lhs.x * rhs.x; }; - friend Number operator*(const int& lhs, const Number& rhs) { Number::nmults++; return lhs * rhs.x; }; - friend Number operator*(const Number& lhs, const int& rhs) { Number::nmults++; return lhs.x * rhs; }; - friend Number operator*(const double& lhs, const Number& rhs) { Number::nmults++; return lhs * rhs.x; }; - friend Number operator*(const Number& lhs, const double& rhs) { Number::nmults++; return lhs.x * rhs; }; - friend Number operator*=(Number& lhs, const Number& rhs) { Number::nmults++; return lhs.x *= rhs.x; }; - friend Number operator/(const Number& lhs, const Number& rhs) { Number::ndivs++; return lhs.x / rhs.x; }; - friend Number operator/(const int& lhs, const Number& rhs) { Number::ndivs++; return lhs / rhs.x; }; - friend Number operator/(const Number& lhs, const int& rhs) { Number::ndivs++; return lhs.x / rhs; }; - friend Number operator/(const double& lhs, const Number& rhs) { Number::ndivs++; return lhs / rhs.x; }; - friend Number operator/(const Number& lhs, const double& rhs) { Number::ndivs++; return lhs.x / rhs; }; - friend Number operator/=(Number& lhs, const Number& rhs) { Number::ndivs++; return lhs.x /= rhs.x; }; - friend bool operator< (const Number& lhs, const Number& rhs){ return lhs.x < rhs.x; } - friend bool operator< (const int& lhs, const Number& rhs){ return lhs < rhs.x; } - friend bool operator< (const Number& lhs, const int& rhs){ return lhs.x < rhs; } - friend bool operator< (const double& lhs, const Number& rhs){ return lhs < rhs.x; } - friend bool operator< (const Number& lhs, const double& rhs){ return lhs.x < rhs; } - friend bool operator> (const Number& lhs, const Number& rhs){ return rhs < lhs; } - friend bool operator> (const int& lhs, const Number& rhs){ return rhs < lhs; } - friend bool operator> (const Number& lhs, const int& rhs){ return rhs < lhs; } - friend bool operator> (const double& lhs, const Number& rhs){ return rhs < lhs; } - friend bool operator> (const Number& lhs, const double& rhs){ return rhs < lhs; } - friend bool operator<=(const Number& lhs, const Number& rhs){ return !(lhs > rhs); } - friend bool operator<=(const int& lhs, const Number& rhs){ return !(lhs > rhs); } - friend bool operator<=(const Number& lhs, const int& rhs){ return !(lhs > rhs); } - friend bool operator<=(const double& lhs, const Number& rhs){ return !(lhs > rhs); } - friend bool operator<=(const Number& lhs, const double& rhs){ return !(lhs > rhs); } - friend bool operator>=(const Number& lhs, const Number& rhs){ return !(lhs < rhs); } - friend bool operator>=(const int& lhs, const Number& rhs){ return !(lhs < rhs); } - friend bool operator>=(const Number& lhs, const int& rhs){ return !(lhs < rhs); } - friend bool operator>=(const double& lhs, const Number& rhs){ return !(lhs < rhs); } - friend bool operator>=(const Number& lhs, const double& rhs){ return !(lhs < rhs); } - friend bool operator==(const Number& lhs, const Number& rhs){ return lhs.x == rhs.x; } - friend bool operator==(const int& lhs, const Number& rhs){ return lhs == rhs.x; } - friend bool operator==(const Number& lhs, const int& rhs){ return lhs.x == rhs; } - friend bool operator==(const double& lhs, const Number& rhs){ return lhs == rhs.x; } - friend bool operator==(const Number& lhs, const double& rhs){ return lhs.x == rhs; } - friend bool operator!=(const Number& lhs, const Number& rhs){ return !(lhs == rhs); } - friend bool operator!=(const int& lhs, const Number& rhs){ return !(lhs == rhs); } - friend bool operator!=(const Number& lhs, const int& rhs){ return !(lhs == rhs); } - friend bool operator!=(const double& lhs, const Number& rhs){ return !(lhs == rhs); } - friend bool operator!=(const Number& lhs, const double& rhs){ return !(lhs == rhs); } - friend std::ostream& operator<<(std::ostream& os, const Number& rhs) { os << rhs.x; return os; } -}; - -#define TSTART -1 -#define TFINAL -2 -#define RESIDUAL -3 -#define ERROR -4 diff --git a/Number.h b/Number.h new file mode 100644 index 0000000..246efa2 --- /dev/null +++ b/Number.h @@ -0,0 +1,31 @@ +#ifndef FPTYPE +#define FPTYPE 2 +#endif + +// 0=half, 1=float, 2=double, 3=quad +#if FPTYPE == 0 +typedef _Float16 fpnumber; +#define FPFMT "%- .4g" +#define FPCAST double +#elif FPTYPE == 1 +typedef float fpnumber; +#define FPFMT "%- .7g" +#define FPCAST double +#elif FPTYPE == 2 +typedef double fpnumber; +#define FPFMT "%- .16g" +#define FPCAST double +#elif FPTYPE == 3 +typedef long double fpnumber; +#define FPFMT ((sizeof(fpnumber)==10)?"%- .19Lg":"%- .34Lg") +#define FPCAST long double +#else +#error UNRECOGNIZED FPTYPE +#endif + +#define Number fpnumber + +#define TSTART -1 +#define TFINAL -2 +#define RESIDUAL -3 +#define ERROR -4 diff --git a/args.C b/args.c similarity index 54% rename from args.C rename to args.c index dc02e1c..e8d5f79 100644 --- a/args.C +++ b/args.c @@ -1,49 +1,76 @@ -#include "heat.H" +#include "heat.h" static char clargs[2048]; -#define HANDLE_ARG(VAR, TYPE, STYLE, HELP) \ +#define HANDLE_HELP(VAR, HELP, STYLE) \ + if (help) \ + {\ + char tmp[256]; \ + int len = snprintf(tmp, sizeof(tmp), " %s=%s(%s%s%s)", \ + VAR, #STYLE, q, valstr, q);\ + snprintf(tmp, sizeof(tmp), "%s", HELP); \ + fprintf(stderr, " %s=%s(%s%s%s)%*s\n", \ + VAR, #STYLE, q, valstr, q, 100-len, tmp);\ + }\ + else \ + { \ + char tmp[64]; \ + fprintf(stderr, " %s=%s(%s%s%s)\n", \ + VAR, #STYLE, q, valstr, q);\ + snprintf(tmp, sizeof(tmp), " %s=%s(%s%s%s)\n", \ + VAR, #STYLE, q, valstr, q);\ + strcat(clargs, tmp); \ + } + +#define HANDLE_FARG(VAR, HELP) \ { \ - char const *style = #STYLE; \ - char const *q = style[1]=='s'?"\"":""; \ void *valp = (void*) &VAR; \ + char const *q = ""; \ int const len = strlen(#VAR)+1; \ - std::stringstream strmvar; \ + char valstr[64]; \ for (i = 1; i < argc; i++) \ {\ - int valid_style = style[1]=='d'||style[1]=='g'||style[1]=='s'; \ if (strncmp(argv[i], #VAR"=", len)) \ continue; \ - assert(valid_style); \ if (strlen(argv[i]+len)) \ - {\ - if (style[1] == 'd') /* int */ \ - *((int*) valp) = (int) strtol(argv[i]+len,0,10); \ - else if (style[1] == 'g') /* double */ \ - *((Number*) valp) = (fpnumber) strtod(argv[i]+len,0); \ - else if (style[1] == 's') /* char* */ \ - *((char**) valp) = (char*) strdup(argv[i]+len); \ - }\ + *((fpnumber*) valp) = (fpnumber) strtod(argv[i]+len,0); \ }\ - strmvar << VAR; \ - if (help) \ + snprintf(valstr, sizeof(valstr), FPFMT, (double) VAR); \ + HANDLE_HELP(#VAR, #HELP, float) \ +} + +#define HANDLE_IARG(VAR, HELP) \ +{ \ + void *valp = (void*) &VAR; \ + char const *q = ""; \ + int const len = strlen(#VAR)+1; \ + char valstr[64]; \ + for (i = 1; i < argc; i++) \ {\ - char tmp[256]; \ - int len = snprintf(tmp, sizeof(tmp), " %s=%s%s%s", \ - #VAR, q, strmvar.str().c_str(), q);\ - snprintf(tmp, sizeof(tmp), "%s (%s)", #HELP, #TYPE); \ - fprintf(stderr, " %s=%s%s%s%*s\n", \ - #VAR, q, strmvar.str().c_str(), q, 80-len, tmp);\ + if (strncmp(argv[i], #VAR"=", len)) \ + continue; \ + if (strlen(argv[i]+len)) \ + *((int*) valp) = (int) strtol(argv[i]+len,0,10); \ }\ - else \ - { \ - char tmp[64]; \ - fprintf(stderr, " %s=%s%s%s\n", \ - #VAR, q, strmvar.str().c_str(), q);\ - snprintf(tmp, sizeof(tmp), " %s=%s%s%s\n", \ - #VAR, q, strmvar.str().c_str(), q);\ - strcat(clargs, tmp); \ - } \ + snprintf(valstr, sizeof(valstr), "%d", (int) VAR); \ + HANDLE_HELP(#VAR, #HELP, int) \ +} + +#define HANDLE_SARG(VAR, HELP) \ +{ \ + void *valp = (void*) &VAR; \ + char const *q = "\""; \ + int const len = strlen(#VAR)+1; \ + char valstr[64]; \ + for (i = 1; i < argc; i++) \ + {\ + if (strncmp(argv[i], #VAR"=", len)) \ + continue; \ + if (strlen(argv[i]+len)) \ + *((char**) valp) = (char*) strdup(argv[i]+len); \ + }\ + snprintf(valstr, sizeof(valstr), "%s", (char*) VAR); \ + HANDLE_HELP(#VAR, #HELP, string) \ } extern Number alpha; @@ -80,11 +107,11 @@ static void handle_ic_help() fprintf(stderr, " ic=\"step(L,Mx,R)\" a step function with value L for x=Mx\n"); fprintf(stderr, " ic=\"rand(S,B,A)\" random values in the range [B-A,B+A] using seed S\n"); fprintf(stderr, " ic=\"sin(A,w)\" a sin wave with amplitude A and frequency w\n"); - fprintf(stderr, " ic=\"spikes(C,A0,X0,A1,X1,...)\" a constant value, C with N spikes of amplitude Ai at position Xi\n"); + fprintf(stderr, " ic=\"spikes(C,A0,X0,A1,X1,...)\" a constant value, C with spikes of amplitude Ai at position Xi\n"); fprintf(stderr, " ic=\"file(foo.dat)\" : read initial condition data from the file foo.dat\n"); fprintf(stderr, "Be sure to use double-quotes (\") as shown and you may also need to set boundary" - "\nconditions such that they *combine* smoothly with the initial conditions.\n"); + "\nconditions such that they *combine* smoothly with the initial condition.\n"); } void @@ -101,27 +128,27 @@ process_args(int argc, char **argv) if (help) fprintf(stderr, "Usage: %s = =...\n", argv[0]); - HANDLE_ARG(alpha, fpnumber, %g, material thermal diffusivity (sq-meters/second)); - HANDLE_ARG(lenx, fpnumber, %g, material length (meters)); - HANDLE_ARG(dx, fpnumber, %g, x-incriment. Best if lenx/dx==int. (meters)); - HANDLE_ARG(dt, fpnumber, %g, t-incriment (seconds)); - HANDLE_ARG(maxt, fpnumber, %g, >0:max sim time (seconds) | <0:min l2 change in soln); - HANDLE_ARG(bc0, fpnumber, %g, boundary condition @ x=0: u(0,t) (Kelvin)); - HANDLE_ARG(bc1, fpnumber, %g, boundary condition @ x=lenx: u(lenx,t) (Kelvin)); - HANDLE_ARG(runame, char*, %s, name to give run and results dir); - HANDLE_ARG(ic, char*, %s, initial condition @ t=0: u(x,0) (Kelvin)); - HANDLE_ARG(alg, char*, %s, algorithm ftcs|dufrank|crankn); + HANDLE_FARG(alpha, material thermal diffusivity (sq-meters/second)); + HANDLE_FARG(lenx, material length (meters)); + HANDLE_FARG(dx, x-incriment. Best if lenx/dx==int. (meters)); + HANDLE_FARG(dt, t-incriment (seconds)); + HANDLE_FARG(maxt, >0:max sim time (seconds) | <0:min l2 change in soln); + HANDLE_FARG(bc0, boundary condition @ x=0: u(0,t) (Kelvin)); + HANDLE_FARG(bc1, boundary condition @ x=lenx: u(lenx,t) (Kelvin)); + HANDLE_SARG(runame, name to give run and results dir); + HANDLE_SARG(ic, initial condition @ t=0: u(x,0) (Kelvin)); + HANDLE_SARG(alg, algorithm ftcs|dufrank|crankn); #ifdef _OPENMP - HANDLE_ARG(nt, int, %d, number of parallel tasks); + HANDLE_IARG(nt, number of parallel tasks); #else - HANDLE_ARG(nt, int, %d, parallel tasking is DISABLED!!); + HANDLE_IARG(nt, parallel tasking is DISABLED!!); nt = 0; #endif - HANDLE_ARG(savi, int, %d, save every i-th solution step); - HANDLE_ARG(save, int, %d, save error in every saved solution); - HANDLE_ARG(outi, int, %d, output progress every i-th solution step); - HANDLE_ARG(noout, int, %d, disable all file outputs); - HANDLE_ARG(prec, const int, %d, precision 0=half/1=float/2=double/3=long double) + HANDLE_IARG(savi, save every i-th solution step); + HANDLE_IARG(save, save error in every saved solution); + HANDLE_IARG(outi, output progress every i-th solution step); + HANDLE_IARG(noout, disable all file outputs); + HANDLE_IARG(prec, precision 1=float/2=double/3=long double) if (help) handle_help(argv[0]); diff --git a/crankn.C b/crankn.c similarity index 95% rename from crankn.C rename to crankn.c index 59809cb..e76e7aa 100644 --- a/crankn.C +++ b/crankn.c @@ -1,4 +1,4 @@ -#include "heat.H" +#include "heat.h" // Licensing: This code is distributed under the GNU LGPL license. // Modified: 30 May 2009 Author: John Burkardt @@ -31,7 +31,7 @@ initialize_crankn(int n, Number const w = alpha * dt / dx / dx; // Build a tri-diagonal matrix - Number *cn_Amat = new Number[3*n](); + Number *cn_Amat = (Number*) malloc(3*n*sizeof(Number)); cn_Amat[0+0*3] = 0.0; cn_Amat[1+0*3] = 1.0; @@ -79,7 +79,7 @@ r83_np_sl ( int n, Number const *a_lu, Number const *b, Number *x) } } -bool +int update_solution_crankn(int n, Number *curr, Number const *last, Number const *cn_Amat, @@ -90,5 +90,5 @@ update_solution_crankn(int n, curr[0] = bc_0; curr[n-1] = bc_1; - return true; + return 1; } diff --git a/dufrank.C b/dufrank.c similarity index 88% rename from dufrank.C rename to dufrank.c index cb2a063..b551fb1 100644 --- a/dufrank.C +++ b/dufrank.c @@ -1,6 +1,6 @@ -#include "heat.H" +#include "heat.h" -bool // false if unstable, true otherwise +int // 0 if unstable, 1 otherwise update_solution_dufrank( int n, // number of samples Number *uk, // new array of u(x,k) to compute/return @@ -22,5 +22,5 @@ update_solution_dufrank( uk[0 ] = bc0; uk[n-1] = bc1; - return true; + return 1; } diff --git a/exact.C b/exact.c similarity index 94% rename from exact.C rename to exact.c index 385cd3a..6610aa9 100644 --- a/exact.C +++ b/exact.c @@ -1,4 +1,4 @@ -#include "heat.H" +#include "heat.h" void compute_exact_steady_state_solution(int n, Number *a, Number dx, char const *ic, diff --git a/ftcs.C b/ftcs.c similarity index 84% rename from ftcs.C rename to ftcs.c index bf20a4e..38c18e9 100644 --- a/ftcs.C +++ b/ftcs.c @@ -1,6 +1,6 @@ -#include "heat.H" +#include "heat.h" -bool // false if unstable, true otherwise +int // false if unstable, true otherwise update_solution_ftcs( int n, // number of samples Number *uk, // new array of u(x,k) to compute/return @@ -12,7 +12,7 @@ update_solution_ftcs( Number r = alpha * dt / (dx * dx); // sanity check for stability - if (r > 0.5) return false; + if (r > 0.5) return 0; // FTCS update algorithm #pragma omp parallel for @@ -23,5 +23,5 @@ update_solution_ftcs( uk[0 ] = bc0; uk[n-1] = bc1; - return true; + return 1; } diff --git a/heat.C b/heat.c similarity index 83% rename from heat.C rename to heat.c index fd3d66a..24a1aee 100644 --- a/heat.C +++ b/heat.c @@ -1,12 +1,6 @@ -#include +#include -#include "heat.H" - -// Number class' statics -int Number::nadds = 0; -int Number::nmults = 0; -int Number::ndivs = 0; -std::size_t Number::nbytes = 0; +#include "heat.h" // Command-line argument variables int noout = 0; @@ -17,14 +11,14 @@ int nt = 0; // number of parallel tasks char const *runame = "heat_results"; char const *alg = "ftcs"; char const *ic = "const(1)"; -Number lenx = Number(1.0); -Number alpha = Number(0.2); -Number dt = Number(0.004); -Number dx = Number(0.1); -Number bc0 = Number(0.0); -Number bc1 = Number(1.0); -Number maxt = Number(2.0); -Number min_change = Number(1e-8*1e-8); +Number lenx = 1.0; +Number alpha = 0.2; +Number dt = 0.004; +Number dx = 0.1; +Number bc0 = 0.0; +Number bc1 = 1.0; +Number maxt = 2.0; +Number min_change = 1e-8*1e-8; // Various arrays of numerical data Number *curr = 0; // current solution @@ -64,19 +58,19 @@ extern void compute_exact_steady_state_solution(int n, Number *a, Number dx, char const *ic, Number alpha, Number t, Number bc0, Number bc1); -extern bool +extern int update_solution_ftcs(int n, Number *curr, Number const *back1, Number alpha, Number dx, Number dt, Number bc_0, Number bc_1); -extern bool +extern int update_solution_crankn(int n, Number *curr, Number const *back1, Number const *cn_Amat, Number bc_0, Number bc_1); -extern bool +extern int update_solution_dufrank(int n, Number *curr, Number const *back1, Number const *back2, Number alpha, Number dx, Number dt, @@ -93,13 +87,13 @@ initialize(void) Nt = (int) (maxt/dt); dx = lenx/(Nx-1); - curr = new Number[Nx](); - back1 = new Number[Nx](); + curr = (Number*) malloc(Nx * sizeof(Number)); + back1 = (Number*) malloc(Nx * sizeof(Number)); if (save) { - exact = new Number[Nx](); - change_history = new Number[Nx](); - error_history = new Number[Nx](); + exact = (Number*) malloc(Nx * sizeof(Number)); + change_history = (Number*) malloc(Nx * sizeof(Number)); + error_history = (Number*) malloc(Nx * sizeof(Number)); } assert(strncmp(alg, "ftcs", 4)==0 || @@ -121,7 +115,7 @@ initialize(void) initialize_crankn(Nx, alpha, dx, dt, &cn_Amat); if (!strncmp(alg, "dufrank", 7)) - back2 = new Number[Nx](); + back2 = (Number*) malloc(Nx * sizeof(Number)); // Initial condition set_initial_condition(Nx, back1, dx, ic); @@ -141,25 +135,22 @@ int finalize(int ti, Number maxt, Number change) if (outi) { printf("Iteration %04d: last change l2=%g\n", ti, (double) change); -#ifndef _OPENMP - printf("Counts: %s\n", Number::counts_string()); -#endif } - delete [] curr; - delete [] back1; - if (back2) delete [] back2; - if (exact) delete [] exact; - if (change_history) delete [] change_history; - if (error_history) delete [] error_history; - if (cn_Amat) delete [] cn_Amat; + free(curr); + free(back1); + if (back2) free(back2); + if (exact) free(exact); + if (change_history) free(change_history); + if (error_history) free(error_history); + if (cn_Amat) free(cn_Amat); if (strncmp(alg, "ftcs", 4)) free((void*)alg); if (strncmp(ic, "const(1)", 8)) free((void*)ic); return retval; } -static bool +static int update_solution() { if (!strcmp(alg, "ftcs")) @@ -168,7 +159,7 @@ update_solution() return update_solution_crankn(Nx, curr, back1, cn_Amat, bc0, bc1); else if (!strcmp(alg, "dufrank")) return update_solution_dufrank(Nx, curr, back1, back2, alpha, dx, dt, bc0, bc1); - return false; + return 0; } static Number diff --git a/heat.H b/heat.h similarity index 61% rename from heat.H rename to heat.h index d0580d7..51b023f 100644 --- a/heat.H +++ b/heat.h @@ -3,19 +3,18 @@ #include #include -#include +#include #ifdef HAVE_FEENABLEEXCEPT #define _GNU_SOURCE #include #endif -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include -#include #ifdef _OPENMP #include @@ -23,5 +22,4 @@ // Lets overlook the fact that we're including one // header inside of another just for convenience -#include "Half.H" -#include "Number.H" +#include "Number.h" diff --git a/helloPy.C b/helloPy.C new file mode 100644 index 0000000..da8d02a --- /dev/null +++ b/helloPy.C @@ -0,0 +1,144 @@ +#include +#include "heat.H" // Include the header for the heat equation solver +// #include // for OpenMP parallelization + +// Define a structure to hold problem data +typedef struct { + double lenx; + double alpha; + int nx; + double dx; + double dt; + double *uk; + double *uk1; +} HeatProblem; + +// Global variable to hold the problem data +static HeatProblem problem; + +/* +// Define a new Python function +static PyObject* foo_func(PyObject *self, PyObject *args) { + return PyUnicode_FromString("hello world"); +} +*/ + +// Function to initialize the heat equation problem +static PyObject* init_problem(PyObject *self, PyObject *args) { + double lenx, alpha; + int nx; + if (!PyArg_ParseTuple(args, "ddi", &lenx, &alpha, &nx)) { + return NULL; + } + + problem.lenx = lenx; + problem.alpha = alpha; + problem.nx = nx; + problem.dx = lenx / (nx - 1); + problem.dt = 0.0; // Default value, will be set in solve function + problem.uk = (double *)malloc(nx * sizeof(double)); + problem.uk1 = (double *)malloc(nx * sizeof(double)); + + if (!problem.uk || !problem.uk1) { + return PyErr_NoMemory(); + } + + // Initialize uk and uk1 with initial conditions (i.e., zero) + for (int i = 0; i < nx; i++) { + problem.uk[i] = 0.0; + problem.uk1[i] = 0.0; + } + + Py_RETURN_NONE; +} + +// Function to solve the heat equation +static PyObject* solve_heat_equation(PyObject *self, PyObject *args) { + double dx, dt, maxt; + int nt; + if (!PyArg_ParseTuple(args, "dddi", &dx, &dt, &maxt, &nt)) { + return NULL; + } + + problem.dx = dx; + problem.dt = dt; + double alpha = problem.alpha; + double lenx = problem.lenx; + int nx = problem.nx; + + double bc0 = 0.0; // Boundary condition at x=0 + double bc1 = 0.0; // Boundary condition at x=lenx + + bool stable = true; + int time_steps = (int)(maxt / dt); + + for (int t = 0; t < time_steps; t++) { + stable = update_solution_ftcs(nx, problem.uk, problem.uk1, alpha, dx, dt, bc0, bc1); + if (!stable) { + PyErr_SetString(PyExc_RuntimeError, "Solution became unstable"); + return NULL; + } + + // Swap uk and uk1 for next iteration + double *temp = problem.uk1; + problem.uk1 = problem.uk; + problem.uk = temp; + } + + // Return the final solution as a Python list + PyObject *result = PyList_New(nx); + for (int i = 0; i < nx; i++) { + PyList_SetItem(result, i, PyFloat_FromDouble(problem.uk[i])); + } + + return result; +} + + +/* +// Declare methods in the module +static PyMethodDef FooMethods[] = { + {"foo_func", foo_func, METH_NOARGS, "Return the string 'hello world'"}, + {NULL, NULL, 0, NULL} // Sentinel +}; +*/ + +// Declare methods in the module +static PyMethodDef HeatMethods[] = { + {"init_problem", init_problem, METH_VARARGS, "Initialize the heat equation problem"}, + {"solve_heat_equation", solve_heat_equation, METH_VARARGS, "Solve the heat equation"}, + {NULL, NULL, 0, NULL} /* Sentinel */ +}; + + +/* +// Define the module structure +static struct PyModuleDef foomodule = { + PyModuleDef_HEAD_INIT, + "foo", // name of module + NULL, // module documentation, may be NULL + -1, // size of per-interpreter state of the module, or -1 if the module keeps state in global variables. + FooMethods +}; +*/ + +// Define the module structure +static struct PyModuleDef heatmodule = { + PyModuleDef_HEAD_INIT, + "heat", /* name of module */ + NULL, /* module documentation, may be NULL */ + -1, /* size of per-interpreter state of the module, or -1 if the module keeps state in global variables. */ + HeatMethods +}; + +/* +// Initialize the module +PyMODINIT_FUNC PyInit_helloPy(void) { + return PyModule_Create(&foomodule); +} +*/ + +// Initialize the module +PyMODINIT_FUNC PyInit_helloPy(void) { + return PyModule_Create(&heatmodule); +} \ No newline at end of file diff --git a/makefile b/makefile index b6633b1..e53f4d9 100644 --- a/makefile +++ b/makefile @@ -11,22 +11,26 @@ RUNAME ?= heat_results PIPEWIDTH ?= 0.1 RM = rm -HDR = Number.H Half.H -SRC = heat.C utils.C args.C exact.C ftcs.C crankn.C dufrank.C -OBJ = $(SRC:.C=.o) -GCOV = $(SRC:.C=.C.gcov) $(SRC:.C=.gcda) $(SRC:.C=.gcno) $(HDR:.H=.H.gcov) +# Headers +HDR = Number.h heat.h +# Source Files +SRC = heat.c utils.c args.c exact.c ftcs.c crankn.c dufrank.c +# Object Files +OBJ = $(SRC:.c=.o) +# Coverage Files +GCOV = $(SRC:.c=.c.gcov) $(SRC:.c=.gcda) $(SRC:.c=.gcno) $(HDR:.h=.h.gcov) +# Executable EXE = heat # Implicit rule for object files -%.o : %.C - $(CXX) -c $(CXXFLAGS) $(CPPFLAGS) $< -o $@ +%.o : %.c + $(CC) -c $(CFLAGS) $(CPPFLAGS) $< -o $@ # Help is default target help: @echo "Targets:" @echo " heat: makes the default heat application (double precision)" @echo " heat-omp: makes default heat application with openmp threads" - @echo " heat-half: makes the heat application with half precision" @echo " heat-single: makes the heat application with single precision" @echo " heat-double: makes the heat application with double precision" @echo " heat-long-double: makes the heat application with long-double precision" @@ -36,10 +40,10 @@ help: # Linking the final heat app heat: $(OBJ) - $(CXX) -o heat $(OBJ) $(LDFLAGS) -lm + $(CC) -o heat $(OBJ) $(LDFLAGS) -lm -heat-omp: CXX=clang -heat-omp: CXXFLAGS=-fopenmp +heat-omp: CC=clang +heat-omp: CFLAGS=-fopenmp heat-omp: LDFLAGS=-lomp -lstdc++ heat-omp: $(OBJ) heat-omp: heat @@ -50,6 +54,7 @@ $(OBJ): $(HDR) # Convenience variable/target for half-precision heat-half: CPPFLAGS=-DFPTYPE=0 +heat-half: CFLAGS=-Wno-format heat-half: $(OBJ) heat-half: heat mv heat heat-half diff --git a/utils.C b/utils.c similarity index 91% rename from utils.C rename to utils.c index 55f1ac8..4bd591b 100644 --- a/utils.C +++ b/utils.c @@ -1,4 +1,4 @@ -#include "heat.H" +#include "heat.h" extern int Nx; extern Number *exact; @@ -77,19 +77,7 @@ write_array(int t, int n, Number dx, Number const *a) outf = fopen(fname,"w"); fprintf(outf, "# %s\n", vname); for (i = 0; i < n; i++) - { -#if FPTYPE == 0 - fprintf(outf, "%- 7.5e %- 7.5e\n", double(i*dx), double(a[i])); -#elif FPTYPE == 1 - fprintf(outf, "%- 11.9e %- 11.9e\n", (double) (i*dx), (double) a[i]); -#elif FPTYPE == 2 - fprintf(outf, "%- 19.17e %- 19.17e\n", (double) (i*dx), (double) a[i]); -#elif FPTYPE == 3 - fprintf(outf, "%- 27.25Le %- 27.25Le\n", (fpnumber) (i*dx), (fpnumber) a[i]); -#elif -#error UNKNOWN FPTYPE -#endif - } + fprintf(outf, FPFMT " " FPFMT "\n", (FPCAST) i*dx, (FPCAST) a[i]); fclose(outf); }