From d220e07aabaf86136e376d73c28a5fe38476a09a Mon Sep 17 00:00:00 2001
From: kar <kas2020@protonmail.com>
Date: Thu, 21 Apr 2022 17:11:53 -0400
Subject: [PATCH] added support for custom values and more modularity
---
src/Science.hpp | 2 +
src/Science/QSPCount.cpp | 55 +++++++++++++++++++++------
src/cases/qsp/qsp.cpp | 82 +++++++++++++++++++++++++++-------------
src/cases/qsp/spins.conf | 47 +++++++++++++++++++++--
4 files changed, 145 insertions(+), 41 deletions(-)
diff --git a/src/Science.hpp b/src/Science.hpp
index f429ffa..cd5c38f 100644
--- a/src/Science.hpp
+++ b/src/Science.hpp
@@ -272,6 +272,8 @@ struct QSPData {
TArrayn::DTArray *temp;
TArrayn::DTArray *rho;
TArrayn::DTArray *salinity;
+ TArrayn::DTArray *custom_T1;
+ TArrayn::DTArray *custom_S1;
TArrayn::DTArray *xgrid;
TArrayn::DTArray *ygrid;
TArrayn::DTArray *zgrid;
diff --git a/src/Science/QSPCount.cpp b/src/Science/QSPCount.cpp
index 358eeaa..25fdae0 100644
--- a/src/Science/QSPCount.cpp
+++ b/src/Science/QSPCount.cpp
@@ -66,6 +66,7 @@ enum QSPType {
QSP_temp,
QSP_rho,
QSP_salinity,
+ QSP_custom,
};
void QSP_write(int local_rank, const QSPVector &local_hist,
@@ -121,6 +122,8 @@ QSPType QSPConvert(const std::string &name) {
converted_type = QSP_rho;
} else if (name.compare("salinity") == 0) {
converted_type = QSP_salinity;
+ } else if (name.compare("custom") == 0) {
+ converted_type = QSP_custom;
}
return converted_type;
}
@@ -149,11 +152,15 @@ TArrayn::DTArray *QSPPtr(const QSPData &qsp_data, const QSPType &type) {
case QSP_salinity:
ptr = qsp_data.salinity;
break;
+ case QSP_custom: // Make sure to set the pointer yourself. Can't do it here.
+ break;
}
return ptr;
}
+// compute the minimum and maximum values for either tracer, then substitute
+// any default values of +- infinity with the global min/max values instead.
void QSPMaxMin(const QSPType &T1_type, const QSPType &S1_type,
TArrayn::DTArray *T1_ptr, TArrayn::DTArray *S1_ptr,
QSPOptions &qsp_options, const QSPData &qsp_data, int i_low,
@@ -214,30 +221,46 @@ void QSPMaxMin(const QSPType &T1_type, const QSPType &S1_type,
MPI_COMM_WORLD);
switch (T1_type) {
case QSP_ke:
- qsp_options.T1_max = glob_ke_max;
- qsp_options.T1_min = glob_ke_min;
+ qsp_options.T1_max =
+ qsp_options.T1_max == double_max ? glob_ke_max : qsp_options.T1_max;
+ qsp_options.T1_min =
+ qsp_options.T1_min == -double_max ? glob_ke_min : qsp_options.T1_min;
break;
case QSP_rho:
- qsp_options.T1_max = glob_rho_max;
- qsp_options.T1_min = glob_rho_min;
+ qsp_options.T1_max =
+ qsp_options.T1_max == double_max ? glob_rho_max : qsp_options.T1_max;
+ qsp_options.T1_min =
+ qsp_options.T1_min == -double_max ? glob_rho_min : qsp_options.T1_min;
break;
default:
- qsp_options.T1_max = psmax(max(*T1_ptr));
- qsp_options.T1_min = psmin(min(*T1_ptr));
+ qsp_options.T1_max = qsp_options.T1_max == double_max
+ ? psmax(max(*T1_ptr))
+ : qsp_options.T1_max;
+ qsp_options.T1_min = qsp_options.T1_min == -double_max
+ ? psmin(min(*T1_ptr))
+ : qsp_options.T1_min;
break;
}
switch (S1_type) {
case QSP_ke:
- qsp_options.S1_max = glob_ke_max;
- qsp_options.S1_min = glob_ke_min;
+ qsp_options.S1_max =
+ qsp_options.S1_max == double_max ? glob_ke_max : qsp_options.S1_max;
+ qsp_options.S1_min =
+ qsp_options.S1_min == -double_max ? glob_ke_min : qsp_options.S1_min;
break;
case QSP_rho:
- qsp_options.S1_max = glob_rho_max;
- qsp_options.S1_min = glob_rho_min;
+ qsp_options.S1_max =
+ qsp_options.S1_max == double_max ? glob_rho_max : qsp_options.S1_max;
+ qsp_options.S1_min =
+ qsp_options.S1_min == -double_max ? glob_rho_min : qsp_options.S1_min;
break;
default:
- qsp_options.S1_max = psmax(max(*S1_ptr));
- qsp_options.S1_min = psmin(min(*S1_ptr));
+ qsp_options.S1_max = qsp_options.S1_max == double_max
+ ? psmax(max(*S1_ptr))
+ : qsp_options.S1_max;
+ qsp_options.S1_min = qsp_options.S1_min == -double_max
+ ? psmin(min(*S1_ptr))
+ : qsp_options.S1_min;
break;
}
} else { // !(cond1 || cond2) == !cond1 && !cond2
@@ -259,6 +282,14 @@ void QSPCount(QSPOptions qsp_options, QSPData qsp_data) {
TArrayn::DTArray *S1_ptr = QSPPtr(qsp_data, S1_type);
TArrayn::DTArray *T1_ptr = QSPPtr(qsp_data, T1_type);
+ if (T1_type == QSP_custom) {
+ T1_ptr = qsp_data.custom_T1;
+ }
+
+ if (S1_type == QSP_custom) {
+ S1_ptr = qsp_data.custom_S1;
+ }
+
if ((!S1_ptr && S1_type != QSP_ke) || (!T1_ptr && T1_type != QSP_ke)) {
std::cout << "Not enough data was provided for the requested tracer. "
"Aborting...\n";
diff --git a/src/cases/qsp/qsp.cpp b/src/cases/qsp/qsp.cpp
index f2cccdf..3dca040 100644
--- a/src/cases/qsp/qsp.cpp
+++ b/src/cases/qsp/qsp.cpp
@@ -8,10 +8,9 @@
#include "../../Options.hpp" // config-file parser
#include "../../Science.hpp" // Science content
-// Defines limits of intrinsic types. Used as default values for
-// T1_max, T1_min, S1_max and/orS1_min
#include <cassert>
#include <limits>
+#include <string>
// Tensor variables for indexing
blitz::firstIndex ii;
@@ -36,21 +35,23 @@ const int z_ind = 2;
// QSP variables
double T1_max, S1_max, T1_min, S1_min;
-string T1_name, S1_name;
+std::string T1_name, S1_name;
int NT, NS;
-string QSP_filename;
-bool use_salinity;
+std::string QSP_filename;
+bool use_salinity, read_rho;
+std::string salinity_filename, temp_filename, rho_filename;
+std::string custom_T1_filename, custom_S1_filename;
+const double double_max = std::numeric_limits<double>::max();
// current output number
int plotnum;
// Derivative options
-string deriv_filenames; // file name to take derivative of
-int start_sequence; // output number to start taking derivatives at
-int final_sequence; // output number to stop taking derivatives at
-int step_sequence; // step between outputs to take derivatives
- // streching/tilting?
-bool v_exist; // Does the v field exist?
+int start_sequence; // output number to start taking derivatives at
+int final_sequence; // output number to stop taking derivatives at
+int step_sequence; // step between outputs to take derivatives
+ // streching/tilting?
+bool v_exist; // Does the v field exist?
/* ------------------ Adjust the class --------------------- */
@@ -152,22 +153,39 @@ public:
qsp_data.temp = NULL;
qsp_data.rho = NULL;
qsp_data.salinity = NULL;
+ qsp_data.custom_T1 = NULL;
+ qsp_data.custom_S1 = NULL;
+
+ DTArray *arr_salinity, *arr_temperature, *arr_rho, *custom_T1, *custom_S1;
- DTArray *arr_salinity, *arr_temperature, *arr_rho;
if (use_salinity) {
arr_salinity = alloc_array(Nx, Ny, Nz);
- init_tracer_restart("s", *arr_salinity);
+ init_tracer_restart(salinity_filename, *arr_salinity);
qsp_data.salinity = arr_salinity;
}
+
+ if (read_rho) {
+ arr_rho = alloc_array(Nx, Ny, Nz);
+ init_tracer_restart(rho_filename, *arr_rho);
+ qsp_data.rho = arr_rho;
+ }
+
if (T1_name.compare("temp") == 0 || S1_name.compare("temp") == 0) {
arr_temperature = alloc_array(Nx, Ny, Nz);
- init_tracer_restart("t", *arr_temperature);
+ init_tracer_restart(temp_filename, *arr_temperature);
qsp_data.temp = arr_temperature;
}
- if (T1_name.compare("rho") == 0 || S1_name.compare("rho") == 0) {
- arr_rho = alloc_array(Nx, Ny, Nz);
- init_tracer_restart("rho", *arr_rho);
- qsp_data.rho = arr_rho;
+
+ if (T1_name.compare("custom") == 0) {
+ custom_T1 = alloc_array(Nx, Ny, Nz);
+ init_tracer_restart(custom_T1_filename, *custom_T1);
+ qsp_data.custom_T1 = custom_T1;
+ }
+
+ if (S1_name.compare("custom") == 0) {
+ custom_S1 = alloc_array(Nx, Ny, Nz);
+ init_tracer_restart(custom_S1_filename, *custom_S1);
+ qsp_data.custom_S1 = custom_S1;
}
QSPCount(qsp_opts, qsp_data);
@@ -216,27 +234,39 @@ int main(int argc, char **argv) {
add_option("type_y", &ygrid_type, "FOURIER", "Grid type in Y");
add_option("type_z", &zgrid_type, "Grid type in Z");
- option_category("Derivative options");
+ option_category("QSP options");
add_option("start_sequence", &start_sequence,
"Sequence number to start taking derivatives at");
add_option("final_sequence", &final_sequence,
"Sequence number to stop taking derivatives at");
add_option("step_sequence", &step_sequence, 1,
"Step between outputs to take derivatives");
+ add_option("salinity_filename", &salinity_filename,
+ "Base Filename of Salinity data (optional).");
+ add_option("temp_filename", &temp_filename,
+ "Base Filename of temperature data (optional).");
+ add_option("T1_filename", &custom_T1_filename,
+ "Base Filename of custom tracer for T1.");
+ add_option("S1_filename", &custom_S1_filename,
+ "Base Filename of custom tracer for S1.");
add_option("T1", &T1_name, "u",
- "Name of tracer 1 for QSP. Valid values are rho,u,v,w,temp or ke");
+ "Name of tracer 1 for QSP. Valid values are rho,u,v,w,temp, "
+ "salinity or ke");
add_option("S1", &S1_name, "w",
- "Name of tracer 2 for QSP. Valid values are rho,u,v,w,temp or ke");
- add_option("T1_max", &T1_max, std::numeric_limits<double>::max(),
+ "Name of tracer 2 for QSP. Valid values are rho,u,v,w,temp, "
+ "salinity or ke");
+ add_option("T1_max", &T1_max, double_max,
"Maximum explicit bin for T1 in QSP.");
- add_option("T1_min", &T1_min, std::numeric_limits<double>::min(),
+ add_option("T1_min", &T1_min, -double_max,
"Minimum explicit bin for T1 in QSP.");
- add_option("S1_max", &S1_max, std::numeric_limits<double>::max(),
+ add_option("S1_max", &S1_max, double_max,
"Maximum explicit bin for S1 in QSP.");
- add_option("S1_min", &S1_min, std::numeric_limits<double>::min(),
+ add_option("S1_min", &S1_min, -double_max,
"Minimum explicit bin for S1 in QSP.");
add_option("salinity", &use_salinity, false,
- "Should salinity be read in from filename s?.");
+ "Should salinity be read in from a file?.");
+ add_option("read_rho", &read_rho, false,
+ "Should rho be read in from a file?.");
add_option("QSP_filename", &QSP_filename, "QSP_default",
"Filename to save data to. Don't include file extension.");
add_option("NS", &NS, 10, "Number of bins for tracer S");
diff --git a/src/cases/qsp/spins.conf b/src/cases/qsp/spins.conf
index b25d2b6..dd65a3b 100644
--- a/src/cases/qsp/spins.conf
+++ b/src/cases/qsp/spins.conf
@@ -15,14 +15,55 @@ type_z = NO_SLIP
mapped_grid = false
v_exist = false
-# Derivative Options
+# If salinity is needed (directly or indirectly) set this to true
+salinity = false
+salinity_filename = s
+
+# If you already pre-computed rho (the equation of state) and saved it to a
+# file, then set this to true. If this is false, and QSP can find the
+# temperature data but not rho, it will use the equation of state function from
+# Science.hpp, and if it also finds salinity it will add that to the equation
+# of state calculation too.
+read_rho = false
+rho_filename = rho
+
+# Sequence Options
start_sequence = 0
final_sequence = 0
step_sequence = 1
-# QSP Parameters
+# Choose what tracers you want to use. T1 and S1 have no distinguishable
+# meaning, just stay consistent with NT and NS. Do note that if you pick ke
+# then QSP will compute it in-place. If you already pre-computed ke and have a
+# saved file for it, use the custom property explained below instead.
+# Tracer options: rho, temp, salinity, u, v, w, ke
+#
+# Note: If you need to use a custom tracer (i.e. not one of the above), set the
+# tracer name to "custom." If you do this however, it is your responsability to
+# ensure that:
+# a) you provide a valid custom filename too and
+# b) all prepocessing is already completed and the data in the file is exactly
+# what you want to use
T1 = temp
S1 = ke
-QSP_filename = QSP
+
+# Tracer base filenames. These are only relevent parameters IFF you set the
+# according tracer name to "custom" above. Otherwise, these are ignored.
+T1_filename = vorticity
+S1_filename = vorticity
+
+# Set the min/max values of the first/last bin in each dimension respectively.
+# These are optional parameters. If you remove these lines from the config QSP
+# will find the global min/max values for each timestep and use that instead.
+# So if you need fixed values for these make sure to include them!
+T1_min = 1
+T1_max = 1
+S1_min = 1
+S1_max = 1
+
+# Number of bins for T1 and S1 chosen above.
NT = 10
NS = 10
+
+# Base filename to save results
+QSP_filename = QSP
--
GitLab