From 4b532e0902a1b9979f4b635bf1059ecc66ac885a Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 18 Jul 2026 20:27:34 -0700 Subject: [PATCH 1/6] chore(tensors_2D_m): uncomment `private` --- src/formal/tensors_2D_m.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index 9349d52..818d439 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -9,7 +9,7 @@ module tensors_2D_m use julienne_m, only : file_t implicit none - !private + private public :: scalar_2D_t public :: vector_2D_t public :: gradient_2D_t From 82e5e610d4b9ed1ccfdc62acb09c770958289aff Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 18 Jul 2026 20:28:31 -0700 Subject: [PATCH 2/6] feat(scalar_2D_t): integer/scalar_2D_t oparator(*) --- src/formal/scalar_2D_s.F90 | 25 +++++++++++++++++++++++++ src/formal/tensors_2D_m.F90 | 23 ++++++++++++++++++++--- 2 files changed, 45 insertions(+), 3 deletions(-) diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index 2eb5076..84ce5d0 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -77,10 +77,35 @@ call_julienne_assert(lhs_x_rhs%consistent()) end procedure + module procedure scalar_2D_postmultiply_integer + + call_julienne_assert(lhs%consistent()) + + lhs_x_rhs = scalar_2D_t( & + tensor_2D_t( & + points = reshape([points_2D_t(lhs%points_(1,1,1,1)%values_ * rhs)], shape = [1,1,1,1]) & + ,cells = lhs%cells_ & + ,x_min = lhs%x_min_ & + ,x_max = lhs%x_max_ & + ,order = lhs%order_ & + ) & + ,gradient_operator_1D_t( & + k = lhs%order_ & + ,dx = (lhs%x_max_ - lhs%x_min_)/lhs%cells_ & + ,cells = lhs%cells_ & + ) ) + + call_julienne_assert(lhs_x_rhs%consistent()) + end procedure + module procedure scalar_2D_premultiply_double lhs_x_rhs = rhs * lhs end procedure + module procedure scalar_2D_premultiply_integer + lhs_x_rhs = rhs * lhs + end procedure + module procedure scalar_2D_plus_scalar call_julienne_assert(rhs%conformable(lhs)) diff --git a/src/formal/tensors_2D_m.F90 b/src/formal/tensors_2D_m.F90 index 818d439..8dfeef2 100644 --- a/src/formal/tensors_2D_m.F90 +++ b/src/formal/tensors_2D_m.F90 @@ -90,7 +90,8 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x contains generic :: assignment(=) => scalar_2D_assign_divergence generic :: operator(.grad.) => scalar_2D_gradient - generic :: operator(*) => scalar_2D_postmultiply_double, scalar_2D_premultiply_double + generic :: operator(*) => scalar_2D_postmultiply_double, scalar_2D_premultiply_double & + ,scalar_2D_postmultiply_integer, scalar_2D_premultiply_integer generic :: operator(+) => scalar_2D_plus_scalar generic :: values => scalar_2D_values generic :: grid => scalar_2D_grid @@ -104,9 +105,9 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x procedure, non_overridable, private :: scalar_2D_values procedure, non_overridable, private :: scalar_2D_grid procedure, non_overridable, private :: scalar_2D_consistent - procedure, non_overridable, private :: scalar_2D_postmultiply_double + procedure, non_overridable, private :: scalar_2D_postmultiply_double, scalar_2D_postmultiply_integer procedure, non_overridable, private :: scalar_2D_plus_scalar - procedure, non_overridable, private, pass(rhs) :: scalar_2D_premultiply_double + procedure, non_overridable, private, pass(rhs) :: scalar_2D_premultiply_double, scalar_2D_premultiply_integer end type interface scalar_2D_t @@ -328,6 +329,14 @@ pure module function scalar_2D_postmultiply_double(lhs, rhs) result(lhs_x_rhs) type(scalar_2D_t) lhs_x_rhs end function + pure module function scalar_2D_postmultiply_integer(lhs, rhs) result(lhs_x_rhs) + !! Result is product of the scalar_2D_t lhs and double-precision (constant) rhs + implicit none + class(scalar_2D_t), intent(in) :: lhs + integer, intent(in) :: rhs + type(scalar_2D_t) lhs_x_rhs + end function + pure module function scalar_2D_premultiply_double(lhs, rhs) result(lhs_x_rhs) !! Result is product of the scalar_2D_t rhs and double-precision (constant) lhs implicit none @@ -336,6 +345,14 @@ pure module function scalar_2D_premultiply_double(lhs, rhs) result(lhs_x_rhs) type(scalar_2D_t) lhs_x_rhs end function + pure module function scalar_2D_premultiply_integer(lhs, rhs) result(lhs_x_rhs) + !! Result is product of the scalar_2D_t rhs and double-precision (constant) lhs + implicit none + class(scalar_2D_t), intent(in) :: rhs + integer, intent(in) :: lhs + type(scalar_2D_t) lhs_x_rhs + end function + pure module function scalar_2D_plus_scalar(lhs, rhs) result(lhs_plus_rhs) !! Result is product of the scalar_2D_t lhs and double-precision (constant) rhs implicit none From 5778b9b69aefb4b7c1fa823d6cc1f1413ea626c7 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 18 Jul 2026 20:49:12 -0700 Subject: [PATCH 3/6] feat(2D-adv-diff): use 4th-order runge-kutta --- example/2D-advection-diffusion.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90 index 1255204..b985757 100644 --- a/example/2D-advection-diffusion.F90 +++ b/example/2D-advection-diffusion.F90 @@ -65,8 +65,14 @@ program advection_diffusion_2D integer step do step = 1, 500 - associate(s_half => s + (dt/2) * d_dt(s, v)) - s = s + dt * d_dt(s_half, v) + associate(k1 => d_dt(s, v)) + associate(k2 => d_dt(s + (dt/2)*k1, v)) + associate(k3 => d_dt(s + (dt/2)*k2, v)) + associate(k4 => d_dt(s + dt*k3, v)) + s = s + (dt/6)*(k1 + 2*k2 + 2*k3 + k4) + end associate + end associate + end associate end associate end do end block advance_time From 9cd76283cbfcd4dd9a18a407303dfd4836352df8 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 18 Jul 2026 21:30:53 -0700 Subject: [PATCH 4/6] feat(2D-adv-diff): increase diffusion This commit increases the diffusion coefficient D in the example 2D advection/diffusion equation solver. Higher diffusion generates more difference betwen the initial- and final-time surface plots of the scalar dependent variable. --- example/2D-advection-diffusion.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90 index b985757..d128eb1 100644 --- a/example/2D-advection-diffusion.F90 +++ b/example/2D-advection-diffusion.F90 @@ -89,7 +89,7 @@ pure function d_dt(s, v) result(ds_dt) type(scalar_2D_t), intent(in) :: s type(vector_2D_t), intent(in) :: v type(scalar_2D_t) ds_dt - double precision, parameter :: D = 0.2D0 + double precision, parameter :: D = 0.5D0 ds_dt = .div. (D * .grad. s) - .div. (v * s) end function From 0e1a44048c2643dd67c91b75942b14efe66d7155 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 18 Jul 2026 22:27:09 -0700 Subject: [PATCH 5/6] feat(gnuplot): plot two 2D scalars side-by-side --- example/scripts/2D-scalar-field.gnuplot | 86 ++++++++++++++++++++----- 1 file changed, 71 insertions(+), 15 deletions(-) diff --git a/example/scripts/2D-scalar-field.gnuplot b/example/scripts/2D-scalar-field.gnuplot index ed8cab0..a8bfccc 100644 --- a/example/scripts/2D-scalar-field.gnuplot +++ b/example/scripts/2D-scalar-field.gnuplot @@ -1,14 +1,25 @@ # ============================================================================ -# 2D-scalar-field.gnuplot -- surface plot CSV +# 2D-scalar-field.gnuplot -- surface plot CSV(s) # Line 1: column labels # Lines 2+: x, y, z data with blank lines between x-slices -# Usage: gnuplot -e "base_name='velocity-potential'" 2D-scalar-field.gnuplot -# Default: base_name='velocity-potential' +# +# Usage (one plot): +# gnuplot -e "base_name='velocity-potential'" 2D-scalar-field.gnuplot +# +# Usage (two plots side by side, sharing the same x/y/z ranges): +# gnuplot -e "base_name='velocity-potential'; base_name2='velocity-potential-2'" \ +# 2D-scalar-field.gnuplot +# +# Defaults: base_name='velocity-potential', base_name2='' (no second plot) # ============================================================================ -if (!exists("base_name")) base_name = "velocity-potential" +if (!exists("base_name")) base_name = "velocity-potential" +if (!exists("base_name2")) base_name2 = "" -datafile = base_name . ".csv" +two_plots = (base_name2 ne "") + +datafile = base_name . ".csv" +datafile2 = base_name2 . ".csv" set datafile separator "," @@ -22,13 +33,40 @@ set datafile separator "," # which is exactly the garbled numeric title ("0.180...E-34(-3.14..., # 3.14...)") you were seeing. Reading the header straight off disk with # the shell sidesteps that entirely. -get_field(n) = system("head -n 1 " . datafile . " | awk -F',' -v n=" . n . " '{v=$n; gsub(/^[ \\t]+|[ \\t]+$/,\"\",v); print v}'") -xlabel = get_field(1) -ylabel = get_field(2) -zlabel = get_field(3) +get_field(file,n) = system("head -n 1 " . file . " | awk -F',' -v n=" . n . " '{v=$n; gsub(/^[ \\t]+|[ \\t]+$/,\"\",v); print v}'") +xlabel = get_field(datafile,1) +ylabel = get_field(datafile,2) +zlabel = get_field(datafile,3) + +# --- 2. Work out a single set of x/y/z ranges shared by both plots -------- +# "stats" only handles one or two "using" columns at a time, so each +# axis is measured separately. Blank lines (the x-slice separators) +# are simply skipped by stats, so they don't interfere. +stats datafile using 1 nooutput +xmin = STATS_min ; xmax = STATS_max +stats datafile using 2 nooutput +ymin = STATS_min ; ymax = STATS_max +stats datafile using 3 nooutput +zmin = STATS_min ; zmax = STATS_max + +if (two_plots) { + stats datafile2 using 1 nooutput + xmin = (STATS_min < xmin) ? STATS_min : xmin + xmax = (STATS_max > xmax) ? STATS_max : xmax + stats datafile2 using 2 nooutput + ymin = (STATS_min < ymin) ? STATS_min : ymin + ymax = (STATS_max > ymax) ? STATS_max : ymax + stats datafile2 using 3 nooutput + zmin = (STATS_min < zmin) ? STATS_min : zmin + zmax = (STATS_max > zmax) ? STATS_max : zmax +} -# --- 2. Plot --- -set title zlabel . "(" . xlabel . ", " . ylabel . ")" +set xrange [xmin:xmax] +set yrange [ymin:ymax] +set zrange [zmin:zmax] +set cbrange [zmin:zmax] + +# --- 3. Common plot styling (applies to both subplots) -------------------- set xlabel xlabel ; set ylabel ylabel set zlabel zlabel offset 3,0 ; set cblabel zlabel set hidden3d @@ -36,12 +74,30 @@ set pm3d depthorder set palette rgbformulae 33,13,10 set ticslevel 0 ; set key off -set terminal gif size 800,600 -set output base_name . ".gif" - # Header is skipped via a shell "tail" pipe rather than "every ::1", so the # blank lines separating x-slices are preserved (pm3d still needs them) and # no per-slice point gets silently dropped the way "every ::1" was doing. -splot "< tail -n +2 " . datafile . "" using 1:2:3 with pm3d title "" +plotcmd(file) = "< tail -n +2 " . file + +# --- 4. Render one or two panels, side by side, with identical axes ------- +if (two_plots) { + set terminal gif size 1500,700 + set output base_name . "_vs_" . base_name2 . ".gif" + set multiplot layout 1,2 + + set title zlabel . "(" . xlabel . ", " . ylabel . ")\n" . base_name + splot plotcmd(datafile) using 1:2:3 with pm3d title "" + + set title zlabel . "(" . xlabel . ", " . ylabel . ")\n" . base_name2 + splot plotcmd(datafile2) using 1:2:3 with pm3d title "" + + unset multiplot +} else { + set terminal gif size 800,600 + set output base_name . ".gif" + + set title zlabel . "(" . xlabel . ", " . ylabel . ")" + splot plotcmd(datafile) using 1:2:3 with pm3d title "" +} set output # flush and close the file From b702af15f62f9857c60702b2b3cdde9e34e6de20 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sun, 19 Jul 2026 12:54:18 -0700 Subject: [PATCH 6/6] feat(gnuplot): 1 script for scalar & vector plots --- .../scripts/2D-scalars-and-vectors.gnuplot | 177 ++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 example/scripts/2D-scalars-and-vectors.gnuplot diff --git a/example/scripts/2D-scalars-and-vectors.gnuplot b/example/scripts/2D-scalars-and-vectors.gnuplot new file mode 100644 index 0000000..245c83f --- /dev/null +++ b/example/scripts/2D-scalars-and-vectors.gnuplot @@ -0,0 +1,177 @@ +# ============================================================================ +# 2D-scalar-field.gnuplot -- one surface/vector-field GIF per input file +# +# Scalar-field files (base_name, base_name2): +# Line 1: column labels +# Lines 2+: x, y, z data with blank lines between x-slices +# +# Vector-field file (vec_name), optional third file: +# Line 1: column labels (x, y, ) +# Lines 2+: x, y, vx, vy (plain rows, one grid point per line) +# +# Every file supplied produces its own GIF (base_name.gif, base_name2.gif, +# vec_name.gif), each a single plot. All plots that are generated together +# in one run share the same x/y/z ranges and the same 3D view angle, so +# they stay visually comparable even though they're separate files -- and +# the vector plot is drawn as a flat plane at the same tilt as the x-y +# plane under the surface plots, not a top-down 2D map. +# +# Usage (one surface): +# gnuplot -e "base_name='velocity-potential'" 2D-scalar-field.gnuplot +# +# Usage (two surfaces, sharing the same x/y/z ranges): +# gnuplot -e "base_name='velocity-potential'; base_name2='velocity-potential-2'" \ +# 2D-scalar-field.gnuplot +# +# Usage (add a horizontal-plane vector-field GIF): +# gnuplot -e "base_name='velocity-potential'; vec_name='velocity'" \ +# 2D-scalar-field.gnuplot +# +# All three together (three GIFs: two surfaces + vectors): +# gnuplot -e "base_name='velocity-potential'; base_name2='velocity-potential-2'; \ +# vec_name='velocity'" 2D-scalar-field.gnuplot +# +# Optional overrides for the vector plot (auto-computed if left unset): +# vec_stride - plot every Nth grid point in each direction (thins arrows) +# vec_scale - multiplier applied to (vx,vy) before drawing each arrow +# +# Defaults: base_name='velocity-potential', base_name2='' (no 2nd surface), +# vec_name='' (no vector plot) +# ============================================================================ + +if (!exists("base_name")) base_name = "velocity-potential" +if (!exists("base_name2")) base_name2 = "" +if (!exists("vec_name")) vec_name = "" +if (!exists("vec_stride")) vec_stride = -1 # -1 = auto +if (!exists("vec_scale")) vec_scale = -1 # -1 = auto + +two_plots = (base_name2 ne "") +have_vec = (vec_name ne "") + +datafile = base_name . ".csv" +datafile2 = base_name2 . ".csv" +vecfile = vec_name . ".csv" + +set datafile separator "," + +# --- 1. Read column headers from line 1 directly via the shell ------------ +get_field(file,n) = system("head -n 1 " . file . " | awk -F',' -v n=" . n . " '{v=$n; gsub(/^[ \\t]+|[ \\t]+$/,\"\",v); print v}'") +xlabel = get_field(datafile,1) +ylabel = get_field(datafile,2) +zlabel = get_field(datafile,3) +if (have_vec) vlabel = get_field(vecfile,3) + +# Build each surface plot's title from the text after the first hyphen in +# its base name, e.g. "scalar-initial" -> "Passive Scalar Initial +# Concentration". toupper() isn't a gnuplot builtin, so the first-letter +# capitalization is done with a short shell/awk round-trip via system(). +cap_first(s) = (strlen(s) > 0 ? system("echo " . s . " | awk '{print toupper(substr($0,1,1)) substr($0,2)}'") : s) +after_hyphen(s) = (strstrt(s,"-") > 0 ? s[strstrt(s,"-")+1:strlen(s)] : s) +scalar_title(s) = "Passive Scalar " . cap_first(after_hyphen(s)) . " Concentration" + +# --- 2. Work out a single set of x/y/z ranges shared by all panels -------- +stats datafile using 1 nooutput +xmin = STATS_min ; xmax = STATS_max +stats datafile using 2 nooutput +ymin = STATS_min ; ymax = STATS_max +stats datafile using 3 nooutput +zmin = STATS_min ; zmax = STATS_max + +if (two_plots) { + stats datafile2 using 1 nooutput + xmin = (STATS_min < xmin) ? STATS_min : xmin + xmax = (STATS_max > xmax) ? STATS_max : xmax + stats datafile2 using 2 nooutput + ymin = (STATS_min < ymin) ? STATS_min : ymin + ymax = (STATS_max > ymax) ? STATS_max : ymax + stats datafile2 using 3 nooutput + zmin = (STATS_min < zmin) ? STATS_min : zmin + zmax = (STATS_max > zmax) ? STATS_max : zmax +} + +set xrange [xmin:xmax] +set yrange [ymin:ymax] +set zrange [zmin:zmax] +set cbrange [zmin:zmax] + +# --- 2b. Work out arrow thinning/scaling for the vector panel ------------- +if (have_vec) { + stats vecfile using 1 nooutput + vxmin = STATS_min ; vxmax = STATS_max + stats vecfile using 2 nooutput + vymin = STATS_min ; vymax = STATS_max + stats vecfile using (sqrt(column(3)**2 + column(4)**2)) nooutput + vmax = STATS_max + + stats vecfile using (abs(column(2)-vymin)<1e-9 ? column(1) : 1/0) nooutput + nx = STATS_records # grid points per row + + target_arrows = 25 + if (vec_stride > 0) { + stride = vec_stride + } else { + stride = floor(nx/target_arrows + 0.5) + if (stride < 1) { stride = 1 } + } + + row(n) = floor(n/nx) + col(n) = n - row(n)*nx + keep(n) = ( (int(row(n)) % int(stride) == 0) && (int(col(n)) % int(stride) == 0) ) ? 1 : 0 + + if (vec_scale > 0) { + vscale = vec_scale + } else { + grid_spacing = (vxmax-vxmin) / (nx/stride) + vscale = (vxmax > 0) ? 0.8*grid_spacing/vmax : 1 + } +} + +# --- 3. Common plot styling (applies to all panels) ------------------------ +if (!exists("view_rotx")) view_rotx = 60 # same 3D view for every panel, so +if (!exists("view_rotz")) view_rotz = 30 # the vector floor lines up with +set view view_rotx,view_rotz # the surfaces' x-y plane below them + +set xlabel xlabel ; set ylabel ylabel +# Removed z-axis label as requested previously +set hidden3d +set pm3d depthorder +set palette rgbformulae 33,13,10 +set ticslevel 0 ; set key off + +# Header is skipped via a shell "tail" pipe rather than "every ::1", so the +# blank lines separating x-slices are preserved (pm3d still needs them) and +# no per-slice point gets silently dropped the way "every ::1" was doing. +plotcmd(file) = "< tail -n +2 " . file + +# --- 4. Render one GIF per input file (each its own plot, own output) ----- +set terminal gif size 800,650 + +set output base_name . ".gif" +set title scalar_title(base_name) +splot plotcmd(datafile) using 1:2:3 with pm3d title "" +set output + +if (two_plots) { + set output base_name2 . ".gif" + set title scalar_title(base_name2) + splot plotcmd(datafile2) using 1:2:3 with pm3d title "" + set output +} + +if (have_vec) { + # --- Dynamic Title Generation --- + cap_vec_name = (strlen(vec_name) > 0 ? system("echo " . vec_name . " | awk '{print toupper(substr($0,1,1)) substr($0,2)}'") : vec_name) + + # 3. Construct the final title: "Velocity Vector Field" + dynamic_title = cap_vec_name . " Vector Field" + + # 4. Determine the output filename with capitalized name + output_filename = cap_vec_name . ".gif" + + set output output_filename + set title dynamic_title + + splot plotcmd(vecfile) using (keep($0) ? column(1) : 1/0):2:(zmin):(column(3)*vscale):(column(4)*vscale):(0) \ + with vectors filled head size 0.08,20 lw 1.3 lc rgb "#1a5fb4" title "" + set output +}