From 18f43a4cb7ec576ad2251a8cce2c1233a6ad75f5 Mon Sep 17 00:00:00 2001 From: "David M. Rogers" Date: Wed, 9 Jun 2021 22:41:23 -0400 Subject: [PATCH 1/6] Removed global statistics, added enum class to write_array. --- Number.H | 109 +++++------------------------------------------------- args.C | 16 ++++---- dufrank.C | 2 +- heat.C | 32 +++------------- heat.H | 1 + make.inc | 8 ++++ makefile | 10 ++--- utils.C | 51 +++++++++++-------------- utils.H | 11 ++++++ 9 files changed, 70 insertions(+), 170 deletions(-) create mode 100644 make.inc create mode 100644 utils.H diff --git a/Number.H b/Number.H index 6f80c8b..2550b36 100644 --- a/Number.H +++ b/Number.H @@ -8,109 +8,20 @@ // 0=half, 1=float, 2=double, 3=quad #if FPTYPE == 0 -typedef half_float::half fpnumber; +typedef half_float::half Number; +template Number castNum(T x) { + return half_float::half_cast(x); +} #elif FPTYPE == 1 -typedef float fpnumber; +typedef float Number; +template Number castNum(T x) { return (Number) x; } #elif FPTYPE == 2 -typedef double fpnumber; +typedef double Number; +template Number castNum(T x) { return (Number) x; } #elif FPTYPE == 3 -typedef long double fpnumber; +typedef long double Number; +template Number castNum(T x) { return (Number) x; } #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]; - sprintf(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; - - // 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/args.C b/args.C index d7eecc3..8e085c9 100644 --- a/args.C +++ b/args.C @@ -20,7 +20,7 @@ static char clargs[2048]; 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); \ + *((Number*) valp) = (Number) strtod(argv[i]+len,0); \ else if (style[1] == 's') /* char* */ \ *((char**) valp) = (char*) strdup(argv[i]+len); \ }\ @@ -100,13 +100,13 @@ 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(alpha, Number, %g, material thermal diffusivity (sq-meters/second)); + HANDLE_ARG(lenx, Number, %g, material length (meters)); + HANDLE_ARG(dx, Number, %g, x-incriment. Best if lenx/dx==int. (meters)); + HANDLE_ARG(dt, Number, %g, t-incriment (seconds)); + HANDLE_ARG(maxt, Number, %g, >0:max sim time (seconds) | <0:min l2 change in soln); + HANDLE_ARG(bc0, Number, %g, boundary condition @ x=0: u(0,t) (Kelvin)); + HANDLE_ARG(bc1, Number, %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); diff --git a/dufrank.C b/dufrank.C index 1fa6ba8..401941e 100644 --- a/dufrank.C +++ b/dufrank.C @@ -11,7 +11,7 @@ update_solution_dufrank( Number bc0, Number bc1) // boundary conditions @ x=0 & x=Lx { Number r = alpha * dt / (dx * dx); - Number q = 1 / (1+r); + Number q = castNum(1.0) / (castNum(1.0)+r); // FTCS update algorithm for (int i = 1; i < n-1; i++) diff --git a/heat.C b/heat.C index bbda2a3..fbf5e45 100644 --- a/heat.C +++ b/heat.C @@ -2,12 +2,6 @@ #include "heat.H" -// Number class' statics -int Number::nadds = 0; -int Number::nmults = 0; -int Number::ndivs = 0; -std::size_t Number::nbytes = 0; - // Command-line argument variables int noout = 0; int savi = 0; @@ -38,19 +32,6 @@ Number *cn_Amat = 0; // A matrix for Crank-Nicholson int Nx; int Nt; -// Utilities -extern Number -l2_norm(int n, Number const *a, Number const *b); - -extern void -copy(int n, Number *dst, Number const *src); - -extern void -write_array(int t, int n, Number dx, Number const *a); - -extern void -set_initial_condition(int n, Number *a, Number dx, char const *ic); - extern void initialize_crankn(int n, Number alpha, Number dx, Number dt, @@ -125,17 +106,16 @@ int finalize(int ti, Number maxt, Number change) { int retval = 0; - write_array(TFINAL, Nx, dx, curr); + write_array(ArrType::TFINAL, Nx, dx, curr); if (save) { - write_array(RESIDUAL, ti, dt, change_history); - write_array(ERROR, ti, dt, error_history); + write_array(ArrType::RESIDUAL, ti, dt, change_history); + write_array(ArrType::ERROR, ti, dt, error_history); } if (outi) { printf("Iteration %04d: last change l2=%g\n", ti, (double) change); - printf("Counts: %s\n", Number::counts_string()); } delete [] curr; @@ -170,13 +150,13 @@ update_output_files(int ti) if (ti>0 && save) { - compute_exact_solution(Nx, exact, dx, ic, alpha, ti*dt, bc0, bc1); + compute_exact_solution(Nx, exact, dx, ic, alpha, castNum(ti*dt), bc0, bc1); if (savi && ti%savi==0) - write_array(ti, Nx, dx, exact); + write_array(ArrType::EXACT, Nx, dx, exact, ti); } if (ti>0 && savi && ti%savi==0) - write_array(ti, Nx, dx, curr); + write_array(ArrType::STEP, Nx, dx, curr, ti); change = l2_norm(Nx, curr, back1); if (save) diff --git a/heat.H b/heat.H index 8f90bd8..a7cae59 100644 --- a/heat.H +++ b/heat.H @@ -20,3 +20,4 @@ // header inside of another just for convenience #include "Half.H" #include "Number.H" +#include "utils.H" diff --git a/make.inc b/make.inc new file mode 100644 index 0000000..79cca14 --- /dev/null +++ b/make.inc @@ -0,0 +1,8 @@ +# ./heat alpha=8.2e-10 lenx=0.25 dx=0.01 dt=100 maxt=5580000 outi=100 savi=1000 bc0=233.15 bc1=294.261 ic="const(294.261)" +ERRBND ?= 1e-6 +PTOOL ?= visit +RUNAME ?= heat_results +PIPEWIDTH ?= 0.1 +RM = rm + +CXXFLAGS = --std=c++14 diff --git a/makefile b/makefile index b3d1eff..f9286dd 100644 --- a/makefile +++ b/makefile @@ -1,9 +1,5 @@ -# ./heat alpha=8.2e-10 lenx=0.25 dx=0.01 dt=100 maxt=5580000 outi=100 savi=1000 bc0=233.15 bc1=294.261 ic="const(294.261)" -ERRBND ?= 1e-6 -PTOOL ?= visit -RUNAME ?= heat_results -PIPEWIDTH ?= 0.1 -RM = rm +# edit environment variables in make.inc +include make.inc HDR = Number.H Half.H SRC = heat.C utils.C args.C exact.C ftcs.C crankn.C dufrank.C @@ -74,7 +70,7 @@ clean: check_clean # # Run for a long time with random initial condition -# and confirm linear stead-state upon termination +# and confirm linear steady-state upon termination # check/check_soln_final.curve: ./heat runame=check outi=0 maxt=10 ic="rand(0,0.2,2)" diff --git a/utils.C b/utils.C index 0c44e42..cd24c77 100644 --- a/utils.C +++ b/utils.C @@ -10,13 +10,13 @@ Number l2_norm(int n, Number const *a, Number const *b) { int i; - Number sum = 0; + Number sum(0.0); for (i = 0; i < n; i++) { Number diff = a[i] - b[i]; sum += diff * diff; } - return sum / n; + return sum / castNum(n); } void @@ -28,7 +28,7 @@ copy(int n, Number *dst, Number const *src) } void -write_array(int t, int n, Number dx, Number const *a) +write_array(ArrType t, int n, Number dx, Number const *a, int time) { int i; char fname[128]; @@ -37,38 +37,31 @@ write_array(int t, int n, Number dx, Number const *a) if (noout) return; - if (t == TSTART) - { + switch(t) { + case ArrType::TSTART: snprintf(fname, sizeof(fname), "%s/%s_soln_00000.curve", runame, runame); snprintf(vname, sizeof(vname), "Temperature"); - } - else if (t == TFINAL) - { + break; + case ArrType::TFINAL: snprintf(fname, sizeof(fname), "%s/%s_soln_final.curve", runame, runame); snprintf(vname, sizeof(vname), "Temperature"); - } - else if (t == RESIDUAL) - { + break; + case ArrType::RESIDUAL: snprintf(fname, sizeof(fname), "%s/%s_change.curve", runame, runame); snprintf(vname, sizeof(vname), "%s/%s_l2_change", runame, runame); - } - else if (t == ERROR) - { + break; + case ArrType::ERROR: snprintf(fname, sizeof(fname), "%s/%s_error.curve", runame, runame); snprintf(vname, sizeof(vname), "%s/%s_l2", runame, runame); - } - else - { - if (a == exact) - { - snprintf(fname, sizeof(fname), "%s/%s_exact_%05d.curve", runame, runame, t); - snprintf(vname, sizeof(vname), "exact_temperature"); - } - else - { - snprintf(fname, sizeof(fname), "%s/%s_soln_%05d.curve", runame, runame, t); - snprintf(vname, sizeof(vname), "Temperature"); - } + break; + case ArrType::EXACT: + snprintf(fname, sizeof(fname), "%s/%s_exact_%05d.curve", runame, runame, time); + snprintf(vname, sizeof(vname), "exact_temperature"); + break; + case ArrType::STEP: + snprintf(fname, sizeof(fname), "%s/%s_soln_%05d.curve", runame, runame, time); + snprintf(vname, sizeof(vname), "Temperature"); + break; } @@ -83,7 +76,7 @@ write_array(int t, int n, Number dx, Number const *a) #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]); + fprintf(outf, "%- 27.25Le %- 27.25Le\n", (Number) (i*dx), (Number) a[i]); #elif #error UNKNOWN FPTYPE #endif @@ -186,5 +179,5 @@ set_initial_condition(int n, Number *a, Number dx, char const *ic) fclose(icfile); free(filename); } - write_array(TSTART, Nx, dx, a); + write_array(ArrType::TSTART, Nx, dx, a); } diff --git a/utils.H b/utils.H new file mode 100644 index 0000000..a1f9431 --- /dev/null +++ b/utils.H @@ -0,0 +1,11 @@ +// Utilities + +enum class ArrType { STEP, EXACT, TSTART, TFINAL, RESIDUAL, ERROR }; + +Number l2_norm(int n, Number const *a, Number const *b); + +void copy(int n, Number *dst, Number const *src); + +void write_array(ArrType t, int n, Number dx, Number const *a, int time=0); + +void set_initial_condition(int n, Number *a, Number dx, char const *ic); From 6e788e51c113e0e336fbbb0ba4ab364429092d41 Mon Sep 17 00:00:00 2001 From: "David M. Rogers" Date: Wed, 9 Jun 2021 23:12:47 -0400 Subject: [PATCH 2/6] Moved command-line arguments into Args struct. --- args.C | 46 +++++++++--------------- args.H | 39 +++++++++++++++++++++ heat.C | 106 +++++++++++++++++++++++--------------------------------- utils.C | 68 +++++++++++++++++------------------- utils.H | 7 ++-- 5 files changed, 135 insertions(+), 131 deletions(-) create mode 100644 args.H diff --git a/args.C b/args.C index 8e085c9..5984281 100644 --- a/args.C +++ b/args.C @@ -6,7 +6,7 @@ static char clargs[2048]; { \ char const *style = #STYLE; \ char const *q = style[1]=='s'?"\"":""; \ - void *valp = (void*) &VAR; \ + void *valp = (void*) &ret.VAR; \ int const len = strlen(#VAR)+1; \ std::stringstream strmvar; \ for (i = 1; i < argc; i++) \ @@ -25,7 +25,7 @@ static char clargs[2048]; *((char**) valp) = (char*) strdup(argv[i]+len); \ }\ }\ - strmvar << VAR; \ + strmvar << ret.VAR; \ if (help) \ {\ char tmp[256]; \ @@ -46,22 +46,7 @@ static char clargs[2048]; } \ } -extern Number alpha; -extern Number lenx; -extern Number dx; -extern Number dt; -extern Number maxt; -extern Number bc0; -extern Number bc1; -extern Number min_change; -extern char const *runame; -extern char const *ic; -extern char const *alg; -extern int savi; -extern int save; -extern int outi; -extern int noout; -int const prec = FPTYPE; +//int const prec = FPTYPE; static void handle_help(char const *argv0) { @@ -86,9 +71,8 @@ static void handle_ic_help() "\nconditions such that they *combine* smoothly with the initial conditions.\n"); } -void -process_args(int argc, char **argv) -{ +Args process_args(int argc, char **argv) { + Args ret{}; int i; int help = 0; @@ -114,12 +98,12 @@ process_args(int argc, char **argv) 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_ARG(prec, const int, %d, precision 0=half/1=float/2=double/3=long double) if (help) handle_help(argv[0]); - if (strcasestr(ic, "help")) + if (strcasestr(ret.ic, "help")) { handle_ic_help(); exit(1); @@ -129,24 +113,26 @@ process_args(int argc, char **argv) exit(1); // Handle possible termination by change threshold criterion - if (maxt < 0) + if (ret.maxt < 0) { - min_change = -maxt * -maxt; - maxt = INT_MAX; + ret.min_change = -ret.maxt * -ret.maxt; + ret.maxt = INT_MAX; } // Handle output results dir creation and save of command-line - if (access(runame, F_OK) == 0) + if (access(ret.runame, F_OK) == 0) { - fprintf(stderr, "An entry \"%s\" already exists\n", runame); + fprintf(stderr, "An entry \"%s\" already exists\n", ret.runame); exit(1); } // Make the output dir and save clargs there too - mkdir(runame, S_IRWXU|S_IRWXG|S_IROTH|S_IXOTH); + mkdir(ret.runame, S_IRWXU|S_IRWXG|S_IROTH|S_IXOTH); char fname[128]; - sprintf(fname, "%s/clargs.out", runame); + sprintf(fname, "%s/clargs.out", ret.runame); FILE *outf = fopen(fname, "w"); fprintf(outf, "%s", clargs); fclose(outf); + + return ret; } diff --git a/args.H b/args.H new file mode 100644 index 0000000..823b38e --- /dev/null +++ b/args.H @@ -0,0 +1,39 @@ +#pragma once + +// Command-line argument variables +struct Args { + int noout; + int savi; + int outi; + int save; + char const *runame; + char const *alg; + char const *ic; + Number lenx; + Number alpha; + Number dt; + Number dx; + Number bc0; + Number bc1; + Number maxt; + Number min_change; + + Args() : + noout ( 0 ), + savi ( 0 ), + outi ( 100 ), + save ( 0 ), + runame ( "heat_results" ), + alg ( "ftcs" ), + ic ( "const(1)" ), + lenx ( Number(1.0) ), + alpha ( Number(0.2) ), + dt ( Number(0.004) ), + dx ( Number(0.1) ), + bc0 ( Number(0.0) ), + bc1 ( Number(1.0) ), + maxt ( Number(2.0) ), + min_change ( Number(1e-8*1e-8) ) {} +}; + +Args process_args(int argc, char **argv); diff --git a/heat.C b/heat.C index fbf5e45..f996e19 100644 --- a/heat.C +++ b/heat.C @@ -2,23 +2,6 @@ #include "heat.H" -// Command-line argument variables -int noout = 0; -int savi = 0; -int outi = 100; -int save = 0; -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); - // Various arrays of numerical data Number *curr = 0; // current solution Number *back1 = 0; // solution back 1 step @@ -37,9 +20,6 @@ initialize_crankn(int n, Number alpha, Number dx, Number dt, Number **_cn_Amat); -extern void -process_args(int argc, char **argv); - extern void compute_exact_solution(int n, Number *a, Number dx, char const *ic, Number alpha, Number t, Number bc0, Number bc1); @@ -69,51 +49,51 @@ update_solution_dufrank(int n, Number *curr, Number bc_0, Number bc_1); static void -initialize(void) +initialize(Args &arg) { - Nx = (int) round((double)(lenx/dx))+1; - Nt = (int) (maxt/dt); - dx = lenx/(Nx-1); + Nx = (int) round((double)(arg.lenx/arg.dx))+1; + Nt = (int) (arg.maxt/arg.dt); + arg.dx = arg.lenx/(Nx-1); curr = new Number[Nx](); back1 = new Number[Nx](); - if (save) + if (arg.save) { exact = new Number[Nx](); change_history = new Number[Nx](); error_history = new Number[Nx](); } - assert(strncmp(alg, "ftcs", 4)==0 || - strncmp(alg, "dufrank", 7)==0 || - strncmp(alg, "crankn", 6)==0); + assert(strncmp(arg.alg, "ftcs", 4)==0 || + strncmp(arg.alg, "dufrank", 7)==0 || + strncmp(arg.alg, "crankn", 6)==0); #ifdef HAVE_FEENABLEEXCEPT feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW); #endif - if (!strncmp(alg, "crankn", 6)) - initialize_crankn(Nx, alpha, dx, dt, &cn_Amat); + if (!strncmp(arg.alg, "crankn", 6)) + initialize_crankn(Nx, arg.alpha, arg.dx, arg.dt, &cn_Amat); - if (!strncmp(alg, "dufrank", 7)) + if (!strncmp(arg.alg, "dufrank", 7)) back2 = new Number[Nx](); /* Initial condition */ - set_initial_condition(Nx, back1, dx, ic); + set_initial_condition(arg, Nx, back1); } -int finalize(int ti, Number maxt, Number change) +int finalize(Args &arg, int ti, Number maxt, Number change) { int retval = 0; - write_array(ArrType::TFINAL, Nx, dx, curr); - if (save) + write_array(ArrType::TFINAL, arg, Nx, arg.dx, curr); + if (arg.save) { - write_array(ArrType::RESIDUAL, ti, dt, change_history); - write_array(ArrType::ERROR, ti, dt, error_history); + write_array(ArrType::RESIDUAL, arg, ti, arg.dt, change_history); + write_array(ArrType::ERROR, arg, ti, arg.dt, error_history); } - if (outi) + if (arg.outi) { printf("Iteration %04d: last change l2=%g\n", ti, (double) change); } @@ -125,41 +105,41 @@ int finalize(int ti, Number maxt, Number change) if (change_history) delete [] change_history; if (error_history) delete [] error_history; if (cn_Amat) delete [] cn_Amat; - if (strncmp(alg, "ftcs", 4)) free((void*)alg); - if (strncmp(ic, "const(1)", 8)) free((void*)ic); + if (strncmp(arg.alg, "ftcs", 4)) free((void*)arg.alg); + if (strncmp(arg.ic, "const(1)", 8)) free((void*)arg.ic); return retval; } static bool -update_solution() +update_solution(Args &arg) { - if (!strcmp(alg, "ftcs")) - return update_solution_ftcs(Nx, curr, back1, alpha, dx, dt, bc0, bc1); - else if (!strcmp(alg, "crankn")) - 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); + if (!strcmp(arg.alg, "ftcs")) + return update_solution_ftcs(Nx, curr, back1, arg.alpha, arg.dx, arg.dt, arg.bc0, arg.bc1); + else if (!strcmp(arg.alg, "crankn")) + return update_solution_crankn(Nx, curr, back1, cn_Amat, arg.bc0, arg.bc1); + else if (!strcmp(arg.alg, "dufrank")) + return update_solution_dufrank(Nx, curr, back1, back2, arg.alpha, arg.dx, arg.dt, arg.bc0, arg.bc1); return false; } static Number -update_output_files(int ti) +update_output_files(Args &arg, int ti) { Number change; - if (ti>0 && save) + if (ti>0 && arg.save) { - compute_exact_solution(Nx, exact, dx, ic, alpha, castNum(ti*dt), bc0, bc1); - if (savi && ti%savi==0) - write_array(ArrType::EXACT, Nx, dx, exact, ti); + compute_exact_solution(Nx, exact, arg.dx, arg.ic, arg.alpha, castNum(ti*arg.dt), arg.bc0, arg.bc1); + if (arg.savi && ti%arg.savi==0) + write_array(ArrType::EXACT, arg, Nx, arg.dx, exact, ti); } - if (ti>0 && savi && ti%savi==0) - write_array(ArrType::STEP, Nx, dx, curr, ti); + if (ti>0 && arg.savi && ti%arg.savi==0) + write_array(ArrType::STEP, arg, Nx, arg.dx, curr, ti); change = l2_norm(Nx, curr, back1); - if (save) + if (arg.save) { change_history[ti] = change; error_history[ti] = l2_norm(Nx, curr, exact); @@ -174,26 +154,26 @@ int main(int argc, char **argv) Number change; // Read command-line args and set values - process_args(argc, argv); + Args arg = process_args(argc, argv); // Allocate arrays and set initial conditions - initialize(); + initialize(arg); // Iterate to max iterations or solution change is below threshold - for (ti = 0; ti*dt < maxt; ti++) + for (ti = 0; ti*arg.dt < arg.maxt; ti++) { // compute the next solution step - if (!update_solution()) + if (!update_solution(arg)) { fprintf(stderr, "Solution criteria violated. Make better choices\n"); exit(1); } // compute amount of change in solution - change = update_output_files(ti); + change = update_output_files(arg, ti); // Handle possible termination by change threshold - if (maxt == INT_MAX && change < min_change) + if (arg.maxt == INT_MAX && change < arg.min_change) { printf("Stopped after %06d iterations for threshold %g\n", ti, (double) change); @@ -201,7 +181,7 @@ int main(int argc, char **argv) } // Output progress - if (outi && ti%outi==0) + if (arg.outi && ti%arg.outi==0) printf("Iteration %04d: last change l2=%g\n", ti, (double) change); // Copy current solution to backi @@ -212,5 +192,5 @@ int main(int argc, char **argv) } // Delete storage and output final results - return finalize(ti, maxt, change); + return finalize(arg, ti, arg.maxt, change); } diff --git a/utils.C b/utils.C index cd24c77..d93eda7 100644 --- a/utils.C +++ b/utils.C @@ -2,8 +2,6 @@ extern int Nx; extern Number *exact; -extern char const *runame; -extern int noout; // Utilities Number @@ -27,39 +25,38 @@ copy(int n, Number *dst, Number const *src) dst[i] = src[i]; } -void -write_array(ArrType t, int n, Number dx, Number const *a, int time) +void write_array(ArrType t, Args &arg, int n, Number step, Number const *a, int time) { int i; char fname[128]; char vname[64]; FILE *outf; - if (noout) return; + if (arg.noout) return; switch(t) { case ArrType::TSTART: - snprintf(fname, sizeof(fname), "%s/%s_soln_00000.curve", runame, runame); + snprintf(fname, sizeof(fname), "%s/%s_soln_00000.curve", arg.runame, arg.runame); snprintf(vname, sizeof(vname), "Temperature"); break; case ArrType::TFINAL: - snprintf(fname, sizeof(fname), "%s/%s_soln_final.curve", runame, runame); + snprintf(fname, sizeof(fname), "%s/%s_soln_final.curve", arg.runame, arg.runame); snprintf(vname, sizeof(vname), "Temperature"); break; case ArrType::RESIDUAL: - snprintf(fname, sizeof(fname), "%s/%s_change.curve", runame, runame); - snprintf(vname, sizeof(vname), "%s/%s_l2_change", runame, runame); + snprintf(fname, sizeof(fname), "%s/%s_change.curve", arg.runame, arg.runame); + snprintf(vname, sizeof(vname), "%s/%s_l2_change", arg.runame, arg.runame); break; case ArrType::ERROR: - snprintf(fname, sizeof(fname), "%s/%s_error.curve", runame, runame); - snprintf(vname, sizeof(vname), "%s/%s_l2", runame, runame); + snprintf(fname, sizeof(fname), "%s/%s_error.curve", arg.runame, arg.runame); + snprintf(vname, sizeof(vname), "%s/%s_l2", arg.runame, arg.runame); break; case ArrType::EXACT: - snprintf(fname, sizeof(fname), "%s/%s_exact_%05d.curve", runame, runame, time); + snprintf(fname, sizeof(fname), "%s/%s_exact_%05d.curve", arg.runame, arg.runame, time); snprintf(vname, sizeof(vname), "exact_temperature"); break; case ArrType::STEP: - snprintf(fname, sizeof(fname), "%s/%s_soln_%05d.curve", runame, runame, time); + snprintf(fname, sizeof(fname), "%s/%s_soln_%05d.curve", arg.runame, arg.runame, time); snprintf(vname, sizeof(vname), "Temperature"); break; } @@ -70,13 +67,13 @@ write_array(ArrType t, int n, Number dx, Number const *a, int time) for (i = 0; i < n; i++) { #if FPTYPE == 0 - fprintf(outf, "%- 7.5e %- 7.5e\n", double(i*dx), double(a[i])); + fprintf(outf, "%- 7.5e %- 7.5e\n", double(i*step), double(a[i])); #elif FPTYPE == 1 - fprintf(outf, "%- 11.9e %- 11.9e\n", (double) (i*dx), (double) a[i]); + fprintf(outf, "%- 11.9e %- 11.9e\n", (double) (i*step), (double) a[i]); #elif FPTYPE == 2 - fprintf(outf, "%- 19.17e %- 19.17e\n", (double) (i*dx), (double) a[i]); + fprintf(outf, "%- 19.17e %- 19.17e\n", (double) (i*step), (double) a[i]); #elif FPTYPE == 3 - fprintf(outf, "%- 27.25Le %- 27.25Le\n", (Number) (i*dx), (Number) a[i]); + fprintf(outf, "%- 27.25Le %- 27.25Le\n", (Number) (i*step), (Number) a[i]); #elif #error UNKNOWN FPTYPE #endif @@ -84,43 +81,42 @@ write_array(ArrType t, int n, Number dx, Number const *a, int time) fclose(outf); } -void -set_initial_condition(int n, Number *a, Number dx, char const *ic) +void set_initial_condition(Args &arg, int n, Number *a) { int i; double x; - if (!strncmp(ic, "const(", 6)) /* const(val) */ + if (!strncmp(arg.ic, "const(", 6)) /* const(val) */ { - double cval = strtod(ic+6, 0); + double cval = strtod(arg.ic+6, 0); for (i = 0; i < n; i++) a[i] = cval; } - else if (!strncmp(ic, "step(", 5)) /* step(left,xmid,right) */ + else if (!strncmp(arg.ic, "step(", 5)) /* step(left,xmid,right) */ { char *p; - double left = strtod(ic+5, &p); + double left = strtod(arg.ic+5, &p); double xmid = strtod(p+1, &p); double right = strtod(p+1, 0); - for (i = 0, x = 0; i < n; i++, x+=dx) + for (i = 0, x = 0; i < n; i++, x+=arg.dx) { if (x < xmid) a[i] = left; else a[i] = right; } } - else if (!strncmp(ic, "ramp(", 5)) /* ramp(left,right) */ + else if (!strncmp(arg.ic, "ramp(", 5)) /* ramp(left,right) */ { char *p; - double left = strtod(ic+5, &p); + double left = strtod(arg.ic+5, &p); double right = strtod(p+1, 0); double dv = (right-left)/(n-1); for (i = 0, x = left; i < n; i++, x+=dv) a[i] = x; } - else if (!strncmp(ic, "rand(", 5)) /* rand(seed,base,amp) */ + else if (!strncmp(arg.ic, "rand(", 5)) /* rand(seed,base,amp) */ { char *p, *ep; - int seed = (int) strtol(ic+5,&p,10); + int seed = (int) strtol(arg.ic+5,&p,10); double base = strtod(p+1, &p); double amp = strtod(p+1, 0); const double maxr = ((long long)1<<31)-1; @@ -128,18 +124,18 @@ set_initial_condition(int n, Number *a, Number dx, char const *ic) for (i = 0; i < n; i++) a[i] = base + amp * (2*random()/maxr - 1); } - else if (!strncmp(ic, "sin(", 4)) /* A*sin(PI*w*x) */ + else if (!strncmp(arg.ic, "sin(", 4)) /* A*sin(PI*w*x) */ { char *p; - double amp = strtod(ic+4,&p); + double amp = strtod(arg.ic+4,&p); double w = strtod(p+1, 0); - for (i = 0, x = 0; i < n; i++, x+=dx) + for (i = 0, x = 0; i < n; i++, x+=arg.dx) a[i] = amp * sin(M_PI*w*x); } - else if (!strncmp(ic, "spikes(", 7)) /* spikes(Const,Amp,Loc,Amp,Loc,...) */ + else if (!strncmp(arg.ic, "spikes(", 7)) /* spikes(Const,Amp,Loc,Amp,Loc,...) */ { char *next; - double cval = strtod(ic+7, &next); + double cval = strtod(arg.ic+7, &next); char const *p = next; for (i = 0, x = 0; i < n; i++) a[i] = cval; @@ -153,10 +149,10 @@ set_initial_condition(int n, Number *a, Number dx, char const *ic) p = ep_idx; } } - else if (!strncmp(ic, "file(", 5)) /* file(numbers.dat) */ + else if (!strncmp(arg.ic, "file(", 5)) /* file(numbers.dat) */ { FILE *icfile; - char *filename = strdup(&ic[5]); + char *filename = strdup(&arg.ic[5]); char *parenchar = strchr(filename, ')'); double val; @@ -179,5 +175,5 @@ set_initial_condition(int n, Number *a, Number dx, char const *ic) fclose(icfile); free(filename); } - write_array(ArrType::TSTART, Nx, dx, a); + write_array(ArrType::TSTART, arg, Nx, arg.dx, a); } diff --git a/utils.H b/utils.H index a1f9431..73be057 100644 --- a/utils.H +++ b/utils.H @@ -1,4 +1,7 @@ // Utilities +#pragma once + +#include "args.H" enum class ArrType { STEP, EXACT, TSTART, TFINAL, RESIDUAL, ERROR }; @@ -6,6 +9,6 @@ Number l2_norm(int n, Number const *a, Number const *b); void copy(int n, Number *dst, Number const *src); -void write_array(ArrType t, int n, Number dx, Number const *a, int time=0); +void write_array(ArrType t, Args &arg, int n, Number dx, Number const *a, int time=0); -void set_initial_condition(int n, Number *a, Number dx, char const *ic); +void set_initial_condition(Args &arg, int n, Number *a); From d3687b50b023ae5e2468a1dace4c66373098d260 Mon Sep 17 00:00:00 2001 From: "David M. Rogers" Date: Thu, 10 Jun 2021 00:05:31 -0400 Subject: [PATCH 3/6] Removed global state. --- Number.H | 24 ++++++++ args.C | 5 ++ args.H | 5 ++ heat.C | 156 ++++++++++++++++++++-------------------------------- heat.H | 1 + solutions.H | 27 +++++++++ utils.C | 18 +++--- utils.H | 6 +- 8 files changed, 136 insertions(+), 106 deletions(-) create mode 100644 solutions.H diff --git a/Number.H b/Number.H index 2550b36..93292a9 100644 --- a/Number.H +++ b/Number.H @@ -1,4 +1,5 @@ #include +#include #include "Half.H" @@ -25,3 +26,26 @@ template Number castNum(T x) { return (Number) x; } #error UNRECOGNIZED FPTYPE #endif +class Vector { + std::vector x; + const int n; + + public: + Vector(int n_) : x(n_), n(n_) { } + + Number& operator[](int i) { + return x[i]; + } + const Number operator[](int i) const { + return x[i]; + } + int size() const { + return n; + } + const Number *data() const { + return x.data(); + } + Number *data() { + return x.data(); + } +}; diff --git a/args.C b/args.C index 5984281..e682c24 100644 --- a/args.C +++ b/args.C @@ -134,5 +134,10 @@ Args process_args(int argc, char **argv) { fprintf(outf, "%s", clargs); fclose(outf); + // final setup on args + ret.Nx = (int) round((double)(ret.lenx/ret.dx))+1; + ret.Nt = (int) (ret.maxt/ret.dt); + ret.dx = ret.lenx/(ret.Nx-1); + return ret; } diff --git a/args.H b/args.H index 823b38e..96672d3 100644 --- a/args.H +++ b/args.H @@ -2,6 +2,10 @@ // Command-line argument variables struct Args { + // Number of points in space, x, and time, t. + int Nx; + int Nt; + int noout; int savi; int outi; @@ -19,6 +23,7 @@ struct Args { Number min_change; Args() : + Nx(0), Nt(0), noout ( 0 ), savi ( 0 ), outi ( 100 ), diff --git a/heat.C b/heat.C index f996e19..510d6a4 100644 --- a/heat.C +++ b/heat.C @@ -2,19 +2,6 @@ #include "heat.H" -// Various arrays of numerical data -Number *curr = 0; // current solution -Number *back1 = 0; // solution back 1 step -Number *back2 = 0; // solution back 2 steps -Number *exact = 0; // exact solution (when available) -Number *change_history = 0; // solution l2norm change history -Number *error_history = 0; // solution error history (when available) -Number *cn_Amat = 0; // A matrix for Crank-Nicholson - -// Number of points in space, x, and time, t. -int Nx; -int Nt; - extern void initialize_crankn(int n, Number alpha, Number dx, Number dt, @@ -24,73 +11,53 @@ extern void compute_exact_solution(int n, Number *a, Number dx, char const *ic, Number alpha, Number t, Number bc0, Number bc1); -extern bool -update_solution_ftcs(int n, - Number *curr, Number const *back1, - Number alpha, Number dx, Number dt, - Number bc_0, Number bc_1); - -extern bool -update_solution_upwind15(int n, - Number *curr, Number const *back1, - Number alpha, Number dx, Number dt, - Number bc_0, Number bc_1); - -extern bool -update_solution_crankn(int n, - Number *curr, Number const *back1, - Number const *cn_Amat, - Number bc_0, Number bc_1); - -extern bool -update_solution_dufrank(int n, Number *curr, - Number const *back1, Number const *back2, - Number alpha, Number dx, Number dt, - Number bc_0, Number bc_1); - -static void -initialize(Args &arg) -{ - Nx = (int) round((double)(arg.lenx/arg.dx))+1; - Nt = (int) (arg.maxt/arg.dt); - arg.dx = arg.lenx/(Nx-1); - - curr = new Number[Nx](); - back1 = new Number[Nx](); - if (arg.save) - { - exact = new Number[Nx](); - change_history = new Number[Nx](); - error_history = new Number[Nx](); - } - - assert(strncmp(arg.alg, "ftcs", 4)==0 || - strncmp(arg.alg, "dufrank", 7)==0 || - strncmp(arg.alg, "crankn", 6)==0); +struct State { + Vector curr ; // current solution + Vector back1 ; // solution back 1 step + Vector back2 ; // solution back 2 steps + Vector exact ; // exact solution (when available) + Vector change_history ; // solution l2norm change history + Vector error_history ; // solution error history (when available) + + Number *cn_Amat ; // A matrix for Crank-Nicholson + + State(Args &arg) + : curr(arg.Nx) + , back1(arg.Nx) + , back2(strncmp(arg.alg, "dufrank", 7)==0 ? arg.Nx : 0) + , exact(arg.save ? arg.Nx : 0) + , change_history((arg.save && arg.Nt != INT_MAX) ? arg.Nt : 0) + , error_history((arg.save && arg.Nt != INT_MAX) ? arg.Nt : 0) + , cn_Amat(nullptr) + { + assert(strncmp(arg.alg, "ftcs", 4)==0 || + strncmp(arg.alg, "dufrank", 7)==0 || + strncmp(arg.alg, "crankn", 6)==0); #ifdef HAVE_FEENABLEEXCEPT - feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW); + feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW); #endif - if (!strncmp(arg.alg, "crankn", 6)) - initialize_crankn(Nx, arg.alpha, arg.dx, arg.dt, &cn_Amat); + if (!strncmp(arg.alg, "crankn", 6)) + initialize_crankn(arg.Nx, arg.alpha, arg.dx, arg.dt, &cn_Amat); - if (!strncmp(arg.alg, "dufrank", 7)) - back2 = new Number[Nx](); - - /* Initial condition */ - set_initial_condition(arg, Nx, back1); -} + /* Initial condition */ + set_initial_condition(arg, back1); + } + ~State() { + if (cn_Amat != nullptr) delete [] cn_Amat; + } +}; -int finalize(Args &arg, int ti, Number maxt, Number change) +int finalize(Args &arg, State &s, int ti, Number maxt, Number change) { int retval = 0; - write_array(ArrType::TFINAL, arg, Nx, arg.dx, curr); + write_array(ArrType::TFINAL, arg, arg.dx, s.curr); if (arg.save) { - write_array(ArrType::RESIDUAL, arg, ti, arg.dt, change_history); - write_array(ArrType::ERROR, arg, ti, arg.dt, error_history); + write_array(ArrType::RESIDUAL, arg, arg.dt, s.change_history, ti); + write_array(ArrType::ERROR, arg, arg.dt, s.error_history, ti); } if (arg.outi) @@ -98,13 +65,7 @@ int finalize(Args &arg, int ti, Number maxt, Number change) printf("Iteration %04d: last change l2=%g\n", ti, (double) change); } - 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; + // FIXME: this free-s an unallocated pointer if the user enters ftcs, etc. if (strncmp(arg.alg, "ftcs", 4)) free((void*)arg.alg); if (strncmp(arg.ic, "const(1)", 8)) free((void*)arg.ic); @@ -112,37 +73,42 @@ int finalize(Args &arg, int ti, Number maxt, Number change) } static bool -update_solution(Args &arg) +update_solution(Args &arg, State &s) { if (!strcmp(arg.alg, "ftcs")) - return update_solution_ftcs(Nx, curr, back1, arg.alpha, arg.dx, arg.dt, arg.bc0, arg.bc1); + return update_solution_ftcs(arg.Nx, s.curr.data(), s.back1.data(), + arg.alpha, arg.dx, arg.dt, arg.bc0, arg.bc1); else if (!strcmp(arg.alg, "crankn")) - return update_solution_crankn(Nx, curr, back1, cn_Amat, arg.bc0, arg.bc1); + return update_solution_crankn(arg.Nx, s.curr.data(), s.back1.data(), + s.cn_Amat, arg.bc0, arg.bc1); else if (!strcmp(arg.alg, "dufrank")) - return update_solution_dufrank(Nx, curr, back1, back2, arg.alpha, arg.dx, arg.dt, arg.bc0, arg.bc1); + return update_solution_dufrank(arg.Nx, s.curr.data(), s.back1.data(), + s.back2.data(), arg.alpha, + arg.dx, arg.dt, arg.bc0, arg.bc1); return false; } static Number -update_output_files(Args &arg, int ti) +update_output_files(Args &arg, State &s, int ti) { Number change; if (ti>0 && arg.save) { - compute_exact_solution(Nx, exact, arg.dx, arg.ic, arg.alpha, castNum(ti*arg.dt), arg.bc0, arg.bc1); + compute_exact_solution(arg.Nx, s.exact.data(), arg.dx, + arg.ic, arg.alpha, castNum(ti*arg.dt), arg.bc0, arg.bc1); if (arg.savi && ti%arg.savi==0) - write_array(ArrType::EXACT, arg, Nx, arg.dx, exact, ti); + write_array(ArrType::EXACT, arg, arg.dx, s.exact, ti); } if (ti>0 && arg.savi && ti%arg.savi==0) - write_array(ArrType::STEP, arg, Nx, arg.dx, curr, ti); + write_array(ArrType::STEP, arg, arg.dx, s.curr, ti); - change = l2_norm(Nx, curr, back1); - if (arg.save) + change = l2_norm(arg.Nx, s.curr.data(), s.back1.data()); + if (arg.save && ti < s.error_history.size()) { - change_history[ti] = change; - error_history[ti] = l2_norm(Nx, curr, exact); + s.change_history[ti] = change; + s.error_history[ti] = l2_norm(arg.Nx, s.curr.data(), s.exact.data()); } return change; @@ -157,20 +123,20 @@ int main(int argc, char **argv) Args arg = process_args(argc, argv); // Allocate arrays and set initial conditions - initialize(arg); + State s(arg); // Iterate to max iterations or solution change is below threshold for (ti = 0; ti*arg.dt < arg.maxt; ti++) { // compute the next solution step - if (!update_solution(arg)) + if (!update_solution(arg, s)) { fprintf(stderr, "Solution criteria violated. Make better choices\n"); exit(1); } // compute amount of change in solution - change = update_output_files(arg, ti); + change = update_output_files(arg, s, ti); // Handle possible termination by change threshold if (arg.maxt == INT_MAX && change < arg.min_change) @@ -185,12 +151,12 @@ int main(int argc, char **argv) printf("Iteration %04d: last change l2=%g\n", ti, (double) change); // Copy current solution to backi - if (back2) - copy(Nx, back2, back1); - copy(Nx, back1, curr); + if (s.back2.size() > 0) + copy(s.back2, s.back1); + copy(s.back1, s.curr); } // Delete storage and output final results - return finalize(arg, ti, arg.maxt, change); + return finalize(arg, s, ti, arg.maxt, change); } diff --git a/heat.H b/heat.H index a7cae59..2b0ab2c 100644 --- a/heat.H +++ b/heat.H @@ -21,3 +21,4 @@ #include "Half.H" #include "Number.H" #include "utils.H" +#include "solutions.H" diff --git a/solutions.H b/solutions.H new file mode 100644 index 0000000..5428c3d --- /dev/null +++ b/solutions.H @@ -0,0 +1,27 @@ +// Declarations for alternative numerical methods to accomplish a time-step. + +#pragma once + +bool +update_solution_ftcs(int n, + Number *curr, Number const *back1, + Number alpha, Number dx, Number dt, + Number bc_0, Number bc_1); + +bool +update_solution_upwind15(int n, + Number *curr, Number const *back1, + Number alpha, Number dx, Number dt, + Number bc_0, Number bc_1); + +bool +update_solution_crankn(int n, + Number *curr, Number const *back1, + Number const *cn_Amat, + Number bc_0, Number bc_1); + +bool +update_solution_dufrank(int n, Number *curr, + Number const *back1, Number const *back2, + Number alpha, Number dx, Number dt, + Number bc_0, Number bc_1); diff --git a/utils.C b/utils.C index d93eda7..e90a2f5 100644 --- a/utils.C +++ b/utils.C @@ -1,8 +1,5 @@ #include "heat.H" -extern int Nx; -extern Number *exact; - // Utilities Number l2_norm(int n, Number const *a, Number const *b) @@ -18,16 +15,18 @@ l2_norm(int n, Number const *a, Number const *b) } void -copy(int n, Number *dst, Number const *src) +copy(Vector &dst, const Vector &src) { + assert(dst.size() == src.size()); int i; - for (i = 0; i < n; i++) + for (i = 0; i < dst.size(); i++) dst[i] = src[i]; } -void write_array(ArrType t, Args &arg, int n, Number step, Number const *a, int time) +void write_array(ArrType t, Args &arg, Number step, const Vector &a, int time) { int i; + int n = a.size(); char fname[128]; char vname[64]; FILE *outf; @@ -46,10 +45,12 @@ void write_array(ArrType t, Args &arg, int n, Number step, Number const *a, int case ArrType::RESIDUAL: snprintf(fname, sizeof(fname), "%s/%s_change.curve", arg.runame, arg.runame); snprintf(vname, sizeof(vname), "%s/%s_l2_change", arg.runame, arg.runame); + n = time < n ? time : n; break; case ArrType::ERROR: snprintf(fname, sizeof(fname), "%s/%s_error.curve", arg.runame, arg.runame); snprintf(vname, sizeof(vname), "%s/%s_l2", arg.runame, arg.runame); + n = time < n ? time : n; break; case ArrType::EXACT: snprintf(fname, sizeof(fname), "%s/%s_exact_%05d.curve", arg.runame, arg.runame, time); @@ -81,8 +82,9 @@ void write_array(ArrType t, Args &arg, int n, Number step, Number const *a, int fclose(outf); } -void set_initial_condition(Args &arg, int n, Number *a) +void set_initial_condition(Args &arg, Vector &a) { + int n = a.size(); int i; double x; @@ -175,5 +177,5 @@ void set_initial_condition(Args &arg, int n, Number *a) fclose(icfile); free(filename); } - write_array(ArrType::TSTART, arg, Nx, arg.dx, a); + write_array(ArrType::TSTART, arg, arg.dx, a); } diff --git a/utils.H b/utils.H index 73be057..5f7668d 100644 --- a/utils.H +++ b/utils.H @@ -7,8 +7,8 @@ enum class ArrType { STEP, EXACT, TSTART, TFINAL, RESIDUAL, ERROR }; Number l2_norm(int n, Number const *a, Number const *b); -void copy(int n, Number *dst, Number const *src); +void copy(Vector &dst, const Vector &src); -void write_array(ArrType t, Args &arg, int n, Number dx, Number const *a, int time=0); +void write_array(ArrType t, Args &arg, Number dx, const Vector &a, int time=0); -void set_initial_condition(Args &arg, int n, Number *a); +void set_initial_condition(Args &arg, Vector &a); From 802e84731348cd99913597b4f1e85fd37c2d5929 Mon Sep 17 00:00:00 2001 From: "David M. Rogers" Date: Thu, 10 Jun 2021 11:21:22 -0400 Subject: [PATCH 4/6] Fixed build flag passthrough. --- make.inc | 2 +- makefile | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/make.inc b/make.inc index 79cca14..4dc83ba 100644 --- a/make.inc +++ b/make.inc @@ -5,4 +5,4 @@ RUNAME ?= heat_results PIPEWIDTH ?= 0.1 RM = rm -CXXFLAGS = --std=c++14 +CSTD = --std=c++14 diff --git a/makefile b/makefile index f9286dd..b880d17 100644 --- a/makefile +++ b/makefile @@ -9,7 +9,7 @@ EXE = heat # Implicit rule for object files %.o : %.C - $(CXX) -c $(CXXFLAGS) $(CPPFLAGS) $< -o $@ + $(CXX) -c $(CSTD) $(CXXFLAGS) $(CPPFLAGS) $< -o $@ # Help is default target help: From 8481750422d2e96c2329f8adf6681beb309557c1 Mon Sep 17 00:00:00 2001 From: "David M. Rogers" Date: Thu, 10 Jun 2021 16:21:43 -0400 Subject: [PATCH 5/6] Changed cn_Amat to Vector type. --- crankn.C | 15 +++++++-------- heat.C | 21 +++------------------ solutions.H | 8 +++++++- 3 files changed, 17 insertions(+), 27 deletions(-) diff --git a/crankn.C b/crankn.C index 59809cb..2046b50 100644 --- a/crankn.C +++ b/crankn.C @@ -4,8 +4,9 @@ // Modified: 30 May 2009 Author: John Burkardt // Modified by Mark C. Miller, July 23, 2017 static void -r83_np_fa(int n, Number *a) +r83_np_fa(Vector &a) { + int n = a.size()/3; int i; for ( i = 1; i <= n-1; i++ ) @@ -22,16 +23,15 @@ r83_np_fa(int n, Number *a) assert( a[1+(n-1)*3] != 0.0 ); } -void +Vector initialize_crankn(int n, - Number alpha, Number dx, Number dt, - Number **_cn_Amat) + Number alpha, Number dx, Number dt) { int i; Number const w = alpha * dt / dx / dx; // Build a tri-diagonal matrix - Number *cn_Amat = new Number[3*n](); + Vector cn_Amat(3*n); cn_Amat[0+0*3] = 0.0; cn_Amat[1+0*3] = 1.0; @@ -49,10 +49,9 @@ initialize_crankn(int n, cn_Amat[2+(n-1)*3] = 0.0; // Factor the matrix. - r83_np_fa(n, cn_Amat); + r83_np_fa(cn_Amat); - // Return the generated matrix - *_cn_Amat = cn_Amat; + return cn_Amat; } // Licensing: This code is distributed under the GNU LGPL license. diff --git a/heat.C b/heat.C index 510d6a4..d4f91ba 100644 --- a/heat.C +++ b/heat.C @@ -2,15 +2,6 @@ #include "heat.H" -extern void -initialize_crankn(int n, - Number alpha, Number dx, Number dt, - Number **_cn_Amat); - -extern void -compute_exact_solution(int n, Number *a, Number dx, char const *ic, - Number alpha, Number t, Number bc0, Number bc1); - struct State { Vector curr ; // current solution Vector back1 ; // solution back 1 step @@ -19,7 +10,7 @@ struct State { Vector change_history ; // solution l2norm change history Vector error_history ; // solution error history (when available) - Number *cn_Amat ; // A matrix for Crank-Nicholson + Vector cn_Amat ; // A matrix for Crank-Nicholson State(Args &arg) : curr(arg.Nx) @@ -28,7 +19,7 @@ struct State { , exact(arg.save ? arg.Nx : 0) , change_history((arg.save && arg.Nt != INT_MAX) ? arg.Nt : 0) , error_history((arg.save && arg.Nt != INT_MAX) ? arg.Nt : 0) - , cn_Amat(nullptr) + , cn_Amat(strncmp(arg.alg, "crankn", 6) == 0 ? initialize_crankn(arg.Nx, arg.alpha, arg.dx, arg.dt) : 0) { assert(strncmp(arg.alg, "ftcs", 4)==0 || strncmp(arg.alg, "dufrank", 7)==0 || @@ -38,15 +29,9 @@ struct State { feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW); #endif - if (!strncmp(arg.alg, "crankn", 6)) - initialize_crankn(arg.Nx, arg.alpha, arg.dx, arg.dt, &cn_Amat); - /* Initial condition */ set_initial_condition(arg, back1); } - ~State() { - if (cn_Amat != nullptr) delete [] cn_Amat; - } }; int finalize(Args &arg, State &s, int ti, Number maxt, Number change) @@ -80,7 +65,7 @@ update_solution(Args &arg, State &s) arg.alpha, arg.dx, arg.dt, arg.bc0, arg.bc1); else if (!strcmp(arg.alg, "crankn")) return update_solution_crankn(arg.Nx, s.curr.data(), s.back1.data(), - s.cn_Amat, arg.bc0, arg.bc1); + s.cn_Amat.data(), arg.bc0, arg.bc1); else if (!strcmp(arg.alg, "dufrank")) return update_solution_dufrank(arg.Nx, s.curr.data(), s.back1.data(), s.back2.data(), arg.alpha, diff --git a/solutions.H b/solutions.H index 5428c3d..8af0285 100644 --- a/solutions.H +++ b/solutions.H @@ -1,7 +1,10 @@ // Declarations for alternative numerical methods to accomplish a time-step. - #pragma once +void +compute_exact_solution(int n, Number *a, Number dx, char const *ic, + Number alpha, Number t, Number bc0, Number bc1); + bool update_solution_ftcs(int n, Number *curr, Number const *back1, @@ -25,3 +28,6 @@ update_solution_dufrank(int n, Number *curr, Number const *back1, Number const *back2, Number alpha, Number dx, Number dt, Number bc_0, Number bc_1); + +Vector +initialize_crankn(int n, Number alpha, Number dx, Number dt); From dc5bf991d62f83ad8a2954ee5d1c2d8f2c2fd07b Mon Sep 17 00:00:00 2001 From: "David M. Rogers" Date: Thu, 10 Jun 2021 16:29:20 -0400 Subject: [PATCH 6/6] Added test for help arg. --- makefile | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/makefile b/makefile index b880d17..96276ae 100644 --- a/makefile +++ b/makefile @@ -98,4 +98,7 @@ check_dufrank: heat check_dufrank/check_dufrank_soln_final.curve cat check_dufrank/check_dufrank_soln_final.curve ./check_lss.sh check_dufrank/check_dufrank_soln_final.curve $(ERRBND) -check_all: check_ftcs check_crankn check_dufrank +check_help: heat + ./heat help 2>&1 | grep -i Examples + +check_all: check_ftcs check_crankn check_dufrank check_help