Commit 59e556b6 authored by Andrew Grace's avatar Andrew Grace Committed by Christopher Subich
Browse files

EOS Scripts

parent bcc8e0ed
......@@ -130,9 +130,126 @@ void R_invt(TArrayn::DTArray & R, TArrayn::DTArray & u,
TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
const string * grid_type, bool v_exist);
inline double nleos_inline(double T, double S){
// Returns the density (kg/m^3)
// MacDougall et. al. 2003 (JAOT 20)
//
// This EOS is a rational function taken at pressure = 0
//
// Constants for Numerator
// First find any regions of negative salinity and set them to zero for this operation.
if (S < 0) {
S = 0;
};
const double N0 = 9.99843699e+02;
const double N1 = 7.35212840e+00;
const double N2 = -5.45928211e-02;
const double N3 = 3.98476704e-04;
const double N4 = 2.96938239e+00;
const double N5 = -7.23268813e-03;
const double N6 = 2.12382341e-03;
const double N7 = 1.04004591e-02;
const double N8 = 1.03970529e-07;
const double N9 = 5.18761880e-06;
const double N10 = -3.24041825e-08;
const double N11 = -1.23869360e-11;
const int p = 0; //Someday we could include pressure effects, and whoever wanted to do that would change p to an input variable from the model
double numerator = N0 + T*(N1 + T*(N2 + N3*T)) + S*(N4 + N5*T + N6*S) + p*(N7 + N8*T*T + N9*S + p*(N10 + N11*T*T));
// Constants for denominator
const double D0 = 1.00000000e+00;
const double D1 = 7.28606739e-03;
const double D2 = -4.60835542e-05;
const double D3 = 3.68390573e-07;
const double D4 = 1.80809186e-10;
const double D5 = 2.14691708e-03;
const double D6 = -9.27062484e-06;
const double D7 = -1.78343643e-10;
const double D8 = 4.76534122e-06;
const double D9 = 1.63410736e-09;
const double D10 = 5.30848875e-06;
const double D11 = -3.03175128e-16;
const double D12 = -1.27934137e-17;
double denominator = D0 + T*(D1 + T*(D2 + T*(D3 + D4*T))) + S*(D5 + T*(D6 + D7*T*T) + sqrt(S)*(D8 + D9*T*T)) + p*(D10 + p*T*(D11*T*T + D12*p));
return numerator/denominator;
}
BZ_DECLARE_FUNCTION2(nleos_inline)
void nleos(TArrayn::DTArray & rho, TArrayn::DTArray & T,
TArrayn::DTArray & S);
inline double compute_alpha(double T0, double S0){
// Computes the thermal expansion coefficient at S0 and T0
// Derivative is a finite difference for simplicity
// Does not divide by rho_0, that is done in lineos()
const double dT = 1e-08;
double TpdT = T0 + dT;
double alpha = (nleos_inline(TpdT,S0) - nleos_inline(T0,S0))/dT;
return alpha;
}
BZ_DECLARE_FUNCTION(compute_alpha)
inline double compute_beta(double T0, double S0){
// Computes the haline contraction coefficient at S0 and T0
// Derivative is a finite difference for simplicity
// Does not divide by rho_0, that is done in lineos()
const double dS = 1e-08;
double SpdS = S0 + dS;
double beta = (nleos_inline(T0,SpdS) - nleos_inline(T0,S0))/dS;
return beta;
}
BZ_DECLARE_FUNCTION(compute_beta)
inline double compute_rho0(double T0, double S0){
// Computes the reference density at S0 and T0
double rho_0 = nleos_inline(T0,S0);
return rho_0;
}
BZ_DECLARE_FUNCTION(compute_rho0)
void lineos(TArrayn::DTArray & rho, TArrayn::DTArray & T,
TArrayn::DTArray & S, const double & T0, const double & S0);
void quadeos(TArrayn::DTArray & rho, TArrayn::DTArray & T);
void eos(const string eos_type, TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S, double T0 = -10, double S0 = -2);
/*
inline double fresh_quad(double T){
// Returns the density (kg/m^3) for water using simple quadratic fit to
// MacDougall et. al. 2003 (JAOT 20)
// Constants are determined via a quadratic fit preserving the Temperature of maximum density
const double rho_max = 9.999744074665388e+02; // Density of freshwater at Tmd (deg C)
//0 pressure relative to sea surface pressure.
const double Tmd = 3.973973973973974; // Temperature of maximum density for freshwater at the surface
const double C = -0.007641729398834; // Fitting Constant
return rho_max + C*(T - Tmd)*(T - Tmd);
}
BZ_DECLARE_FUNCTION(fresh_quad)
*/
// Equation of state for seawater, polynomial fit from
// Brydon, Sun, Bleck (1999) (JGR)
//
inline double eqn_of_state(double T, double S){
// Returns the density anomaly (kg/m^3) for water at the given
// temperature T (degrees celsius) and salinity S (PSU?)
......
#include "../Science.hpp"
#include "math.h"
#include "../Par_util.hpp"
#include "stdio.h"
#include "../Split_reader.hpp"
#include "../T_util.hpp"
#include "../Parformer.hpp"
#include "../Sorter.hpp"
#include <numeric>
#include <stdlib.h>
using blitz::Array;
using blitz::cos;
using namespace TArrayn;
using namespace NSIntegrator;
using namespace Transformer;
void eos(const string eos_type, TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S, double T0, double S0) {
bool is_T0_valid = false;
bool is_S0_valid = false;
if (T0 >= -2) {
is_T0_valid = true;
}
if (S0 >= 0) {
is_S0_valid = true;
}
if (eos_type == "QUADEOS") {
//For Debugging purposes
/*
if (master()) {
fprintf(stdout,"You picked the QUADEOS option.\n");
}
MPI_Finalize(); exit(0);
*/
//Call quadeos
quadeos(rho,T);
}
else if (eos_type == "LINEOS") {
if (is_T0_valid & is_S0_valid) {
//For Debugging purposes
/*
if (master()) {
fprintf(stdout,"You picked the LINEOS option.\n");
}
MPI_Finalize(); exit(0);
*/
//Call lineos
lineos(rho,T,S,T0,S0);
}
else {
if (master()) {
fprintf(stderr,"Invalid option for background temperature or salinity. Exiting.\n");
}
MPI_Finalize(); exit(0);
}
}
else if (eos_type == "NLEOS") {
//For Debugging purposes
/*
if (master()) {
fprintf(stdout,"You picked the NLEOS option.\n");
}
MPI_Finalize(); exit(0);
*/
//Call nleos
nleos(rho,T,S);
}
else {
if (master()) {
fprintf(stderr,"Invalid option received for eos type. Exiting.\n");
}
MPI_Finalize(); exit(0);
}
}
#include "../Science.hpp"
#include "math.h"
#include "../Par_util.hpp"
#include "stdio.h"
#include "../Split_reader.hpp"
#include "../T_util.hpp"
#include "../Parformer.hpp"
#include "../Sorter.hpp"
#include <numeric>
using blitz::Array;
using blitz::cos;
using namespace TArrayn;
using namespace NSIntegrator;
using namespace Transformer;
void lineos(TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S,
const double & T0, const double & S0) {
double alpha = compute_alpha(T0,S0);
double beta = compute_beta(T0,S0);
double rho_0 = compute_rho0(T0,S0);
rho = rho_0*(1 + alpha/rho_0*(T - T0) + beta/rho_0*(S - S0));
}
#include "../Science.hpp"
#include "math.h"
#include "../Par_util.hpp"
#include "stdio.h"
#include "../Split_reader.hpp"
#include "../T_util.hpp"
#include "../Parformer.hpp"
#include "../Sorter.hpp"
#include <numeric>
using blitz::Array;
using blitz::cos;
using namespace TArrayn;
using namespace NSIntegrator;
using namespace Transformer;
void nleos(TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S) {
// Returns the density (kg/m^3)
// MacDougall et. al. 2003 (JAOT 20)
//
// This EOS is a rational function taken at pressure = 0
//
// First find any regions of negative salinity and set them to zero for this operation.
//Numerator
const double N0 = 9.99843699e+02;
const double N1 = 7.35212840e+00;
const double N2 = -5.45928211e-02;
const double N3 = 3.98476704e-04;
const double N4 = 2.96938239e+00;
const double N5 = -7.23268813e-03;
const double N6 = 2.12382341e-03;
const double N7 = 1.04004591e-02;
const double N8 = 1.03970529e-07;
const double N9 = 5.18761880e-06;
const double N10 = -3.24041825e-08;
const double N11 = -1.23869360e-11;
const int p = 0; //Someday we could include pressure effects, and whoever wanted to do that would change p to an input variable from the model
rho = N0 + T*(N1 + T*(N2 + N3*T)) + S*(N4 + N5*T + N6*S) + p*(N7 + N8*T*T + N9*S + p*(N10 + N11*T*T));
// Constants for denominator
const double D0 = 1.00000000e+00;
const double D1 = 7.28606739e-03;
const double D2 = -4.60835542e-05;
const double D3 = 3.68390573e-07;
const double D4 = 1.80809186e-10;
const double D5 = 2.14691708e-03;
const double D6 = -9.27062484e-06;
const double D7 = -1.78343643e-10;
const double D8 = 4.76534122e-06;
const double D9 = 1.63410736e-09;
const double D10 = 5.30848875e-06;
const double D11 = -3.03175128e-16;
const double D12 = -1.27934137e-17;
rho = rho/(D0 + T*(D1 + T*(D2 + T*(D3 + D4*T))) + S*(D5 + T*(D6 + D7*T*T) + sqrt(max(S,0))*(D8 + D9*T*T)) + p*(D10 + p*T*(D11*T*T + D12*p)));
}
#include "../Science.hpp"
#include "math.h"
#include "../Par_util.hpp"
#include "stdio.h"
#include "../Split_reader.hpp"
#include "../T_util.hpp"
#include "../Parformer.hpp"
#include "../Sorter.hpp"
#include <numeric>
using blitz::Array;
using blitz::cos;
using namespace TArrayn;
using namespace NSIntegrator;
using namespace Transformer;
void quadeos(TArrayn::DTArray & rho, TArrayn::DTArray & T) {
// Returns the density (kg/m^3) for water using simple quadratic fit to
// MacDougall et. al. 2003 (JAOT 20)
// Constants are determined via a quadratic fit preserving the Temperature of maximum density
const double rho_max = 9.999744074665388e+02; // Density of freshwater at Tmd (deg C)
//0 pressure relative to sea surface pressure.
const double Tmd = 3.973973973973974; // Temperature of maximum density for freshwater at the surface
const double C = -0.007641729398834; // Fitting Constant
rho = rho_max + C*(T - Tmd)*(T - Tmd);
}
This diff is collapsed.
## Salt-temperature configuration file
# Spatial Parameters
Lx = 1.0
Ly = 1.0
Lz = 0.1
Nx = 1024
Ny = 1
Nz = 128
min_x = 0
min_y = 0
min_z = 0
# Expansion types
type_x = FREE_SLIP
type_y = FREE_SLIP
type_z = FREE_SLIP
mapped_grid = false
# Physical Parameters
g = 9.81
rot_f = 0.0e-3
rho_0 = 1000.0
visco = 2e-6
kappa_T = 2e-6
kappa_S = 1e-7
# Problem Parameters (delta_rho is percentage)
delta_T = 25.0
delta_S = 5.0
T_0 = 15.0
S_0 = 0.0
delta_x = 0.02
Lmix = 0.1
# Problem topography parameters
hill_height = 0.1
hill_centre = 0.5
hill_width = 0.05
# Temporal Parameters
final_time = 100
plot_interval = 1
#dt_max = 0.0
# Restart Options
restart = false
restart_time = 0.0
restart_sequence = 0
restart_from_dump = false
compute_time = -1
# Perturbation Parameter
perturb = 1e-3
# Filter Parameters
f_cutoff = 0.6
f_order = 2.0
f_strength = 20.0
# secondary diagnostics
#compute_enstrophy = true
#compute_dissipation = true
compute_BPE = false
compute_internal_to_BPE = false
#compute_stresses_top = false
#compute_stresses_bottom = false
This diff is collapsed.
## Salt-temperature configuration file
# Spatial Parameters
Lx = 1.0
Ly = 1.0
Lz = 0.1
Nx = 1024
Ny = 1
Nz = 128
min_x = 0
min_y = 0
min_z = 0
# Expansion types
type_x = FREE_SLIP
type_y = FREE_SLIP
type_z = FREE_SLIP
mapped_grid = false
# Physical Parameters
g = 9.81
rot_f = 0.0e-3
rho_0 = 1000.0
visco = 2e-6
kappa_T = 2e-6
kappa_S = 1e-7
# Problem Parameters (delta_rho is percentage)
delta_T = 25.0
delta_S = 5.0
T_0 = 15.0
S_0 = 0.0
delta_x = 0.02
Lmix = 0.1
# Problem topography parameters
hill_height = 0.1
hill_centre = 0.5
hill_width = 0.05
# Temporal Parameters
final_time = 100
plot_interval = 1
#dt_max = 0.0
# Restart Options
restart = false
restart_time = 0.0
restart_sequence = 0
restart_from_dump = false
compute_time = -1
# Perturbation Parameter
perturb = 1e-3
# Filter Parameters
f_cutoff = 0.6
f_order = 2.0
f_strength = 20.0
# secondary diagnostics
#compute_enstrophy = true
#compute_dissipation = true
compute_BPE = false
compute_internal_to_BPE = false
#compute_stresses_top = false
#compute_stresses_bottom = false
This diff is collapsed.
## Salt-temperature configuration file
# Spatial Parameters
Lx = 1.0
Ly = 1.0
Lz = 0.1
Nx = 1024
Ny = 1
Nz = 128
min_x = 0
min_y = 0
min_z = 0
# Expansion types
type_x = FREE_SLIP
type_y = FREE_SLIP
type_z = FREE_SLIP
mapped_grid = false
# Physical Parameters
g = 9.81
rot_f = 0.0e-3
rho_0 = 1000.0
visco = 2e-6
kappa_T = 2e-6
kappa_S = 1e-7
# Problem Parameters (delta_rho is percentage)
delta_T = 25.0
delta_S = 5.0
T_0 = 15.0
S_0 = 0.0
delta_x = 0.02
Lmix = 0.1
# Problem topography parameters
hill_height = 0.1
hill_centre = 0.5
hill_width = 0.05
# Temporal Parameters
final_time = 100
plot_interval = 1
#dt_max = 0.0
# Restart Options
restart = false
restart_time = 0.0
restart_sequence = 0
restart_from_dump = false
compute_time = -1
# Perturbation Parameter
perturb = 1e-3
# Filter Parameters
f_cutoff = 0.6
f_order = 2.0
f_strength = 20.0
# secondary diagnostics
#compute_enstrophy = true
#compute_dissipation = true
compute_BPE = false
compute_internal_to_BPE = false
#compute_stresses_top = false
#compute_stresses_bottom = false
This diff is collapsed.
## Salt-temperature configuration file
# Spatial Parameters
Lx = 1.0
Ly = 1.0
Lz = 0.1
Nx = 1024
Ny = 1
Nz = 128
min_x = 0
min_y = 0
min_z = 0
# Expansion types
type_x = FREE_SLIP
type_y = FREE_SLIP
type_z = FREE_SLIP
mapped_grid = false
# Physical Parameters
g = 9.81
rot_f = 0.0e-3
rho_0 = 1000.0
visco = 2e-6
kappa_T = 2e-6
kappa_S = 1e-7
# Problem Parameters (delta_rho is percentage)
delta_T = 25.0
delta_S = 5.0
T_0 = 15.0
S_0 = 0.0
delta_x = 0.02
Lmix = 0.1
# Problem topography parameters
hill_height = 0.1
hill_centre = 0.5
hill_width = 0.05
# Temporal Parameters
final_time = 100
plot_interval = 1
#dt_max = 0.0
# Restart Options
restart = false
restart_time = 0.0
restart_sequence = 0
restart_from_dump = false
compute_time = -1
# Perturbation Parameter
perturb = 1e-3
# Filter Parameters
f_cutoff = 0.6
f_order = 2.0
f_strength = 20.0
# secondary diagnostics
#compute_enstrophy = true
#compute_dissipation = true
compute_BPE = false
compute_internal_to_BPE = false
#compute_stresses_top = false
#compute_stresses_bottom = false
This diff is collapsed.
This diff is collapsed.