From da66f8eeba30e48620f46877fcc994da33043b3a Mon Sep 17 00:00:00 2001 From: cmkaul Date: Thu, 3 Sep 2015 16:29:31 +0200 Subject: [PATCH] added to output of TKE SGS scheme --- Csrc/kinematics.h | 7 ++++--- Csrc/sgs.h | 7 +++---- Kinematics.pxd | 3 ++- Kinematics.pyx | 9 +++++---- SGS.pyx | 14 ++++++++++++-- Simulation3d.pyx | 4 ++-- 6 files changed, 28 insertions(+), 16 deletions(-) diff --git a/Csrc/kinematics.h b/Csrc/kinematics.h index cc46befe..1a75c891 100644 --- a/Csrc/kinematics.h +++ b/Csrc/kinematics.h @@ -141,7 +141,8 @@ void compute_strain_rate_mag(const struct DimStruct *dims, double* restrict stra return; } -void compute_wind_speed_angle(const struct DimStruct *dims, double* restrict u, double* restrict v, double* restrict wind_speed, double* restrict wind_angle){ +void compute_wind_speed_angle(const struct DimStruct *dims, double* restrict u, double* restrict v, + double* restrict wind_speed, double* restrict wind_angle, double u0, double v0){ const size_t istride = dims->nlg[1] * dims->nlg[2]; const size_t jstride = dims->nlg[2]; @@ -160,8 +161,8 @@ void compute_wind_speed_angle(const struct DimStruct *dims, double* restrict u, const size_t jshift = j*jstride; for(size_t k=kmin;kdx[0]*dims->dx[1]*dims->dx[2]); for (size_t i=0; inpg; i++){ - double s2 = strain_rate_mag[i]*strain_rate_mag[i]; - visc[i] = cs*cs*delta*delta*s2; + visc[i] = cs*cs*delta*delta*strain_rate_mag[i]; if(buoy_freq[i] > 0.0){ - double fb = sqrt(fmax(1.0 - buoy_freq[i]/(prt*s2),0.0)); + double fb = sqrt(fmax(1.0 - buoy_freq[i]/(prt*strain_rate_mag[i]*strain_rate_mag[i]),0.0)); visc[i] = visc[i] * fb; } diff[i] = visc[i]/prt; @@ -21,7 +20,7 @@ double* restrict buoy_freq, double* restrict strain_rate_mag, double cs, double double tke_ell(double cn, double e, double buoy_freq, double delta){ double ell; if(buoy_freq> 1.0e-10){ - ell = fmax(fmin(cn*sqrt(fmax(e,0.0)/buoy_freq),delta),1.0e-10); + ell = fmax(fmin(cn*sqrt(fmax(e,0.0)/buoy_freq),delta),1e-10); } else{ ell = delta; diff --git a/Kinematics.pxd b/Kinematics.pxd index b56b083d..b370e69d 100644 --- a/Kinematics.pxd +++ b/Kinematics.pxd @@ -1,6 +1,7 @@ cimport Grid cimport PrognosticVariables cimport ParallelMPI +cimport ReferenceState from NetCDFIO cimport NetCDFIO_Stats cdef class Kinematics: @@ -13,4 +14,4 @@ cdef class Kinematics: Py_ssize_t get_grad_shift(self, Grid.Grid Gr, Py_ssize_t vel_i, Py_ssize_t dx_j) cpdef initialize(self, Grid.Grid Gr, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa) cpdef update(self, Grid.Grid Gr, PrognosticVariables.PrognosticVariables PV) - cpdef stats_io(self, Grid.Grid Gr, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa) + cpdef stats_io(self, Grid.Grid Gr,ReferenceState.ReferenceState Ref, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa) diff --git a/Kinematics.pyx b/Kinematics.pyx index 6e460096..853183e5 100644 --- a/Kinematics.pyx +++ b/Kinematics.pyx @@ -1,6 +1,7 @@ cimport Grid cimport PrognosticVariables cimport ParallelMPI +cimport ReferenceState from NetCDFIO cimport NetCDFIO_Stats import numpy as np @@ -10,7 +11,7 @@ cdef extern from "kinematics.h": void compute_velocity_gradient(Grid.DimStruct *dims, double *v, double *vgrad, long d) void compute_strain_rate(Grid.DimStruct *dims, double *vgrad, double *strain_rate) void compute_strain_rate_mag(Grid.DimStruct *dims, double *vgrad, double *strain_rate) - void compute_wind_speed_angle(Grid.DimStruct *dims, double *u, double *v, double *wind_speed, double *wind_angle) + void compute_wind_speed_angle(Grid.DimStruct *dims, double *u, double *v, double *wind_speed, double *wind_angle, double u0, double v0) cdef class Kinematics: def __init__(self): @@ -50,7 +51,7 @@ cdef class Kinematics: return - cpdef stats_io(self, Grid.Grid Gr, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa): + cpdef stats_io(self, Grid.Grid Gr, ReferenceState.ReferenceState RS, PrognosticVariables.PrognosticVariables PV, NetCDFIO_Stats NS, ParallelMPI.ParallelMPI Pa): cdef: double [:] mean = np.empty((Gr.dims.nlg[2],),dtype=np.double,order='c') Py_ssize_t u_shift = PV.get_varshift(Gr, 'u') @@ -60,12 +61,12 @@ cdef class Kinematics: mean = Pa.HorizontalMean(Gr,&self.strain_rate_mag[0]) NS.write_profile('strain_rate_magnitude',mean[Gr.dims.gw:-Gr.dims.gw],Pa) - compute_wind_speed_angle(&Gr.dims, &PV.values[u_shift], &PV.values[v_shift], &self.wind_speed[0], &self.wind_angle[0]) + compute_wind_speed_angle(&Gr.dims, &PV.values[u_shift], &PV.values[v_shift], &self.wind_speed[0], &self.wind_angle[0],RS.u0, RS.v0) mean = Pa.HorizontalMean(Gr,&self.wind_speed[0]) NS.write_profile('wind_speed',mean[Gr.dims.gw:-Gr.dims.gw],Pa) - mean = Pa.HorizontalMean(Gr,&self.wind_speed[0]) + mean = Pa.HorizontalMean(Gr,&self.wind_angle[0]) NS.write_profile('wind_angle',mean[Gr.dims.gw:-Gr.dims.gw],Pa) return diff --git a/SGS.pyx b/SGS.pyx index 2aa24f09..a69b9088 100644 --- a/SGS.pyx +++ b/SGS.pyx @@ -17,6 +17,7 @@ cdef extern from "sgs.h": void tke_shear_production(Grid.DimStruct *dims, double* e_tendency, double* visc, double* strain_rate_mag) void tke_buoyant_production(Grid.DimStruct *dims, double* e_tendency, double* diff, double* buoy_freq) void tke_surface(Grid.DimStruct *dims, double* e, double* lmo, double* ustar, double h_bl, double zb) nogil + double tke_ell(double cn, double e, double buoy_freq, double delta) nogil cdef class SGS: def __init__(self,namelist): @@ -151,6 +152,7 @@ cdef class TKE: NS.add_profile('tke_shear_tendency', Gr, Pa) NS.add_profile('tke_buoyancy_tendency', Gr, Pa) NS.add_profile('tke_prandtl_number', Gr, Pa) + NS.add_profile('tke_mixing_length', Gr, Pa) return @@ -221,7 +223,7 @@ cdef class TKE: Py_ssize_t e_shift = PV.get_varshift(Gr,'e') double [:] tmp_tendency = np.zeros((Gr.dims.npg),dtype=np.double,order='c') double [:] mean_tendency = np.empty((Gr.dims.nlg[2],),dtype=np.double,order='c') - double [:] prt = np.zeros((Gr.dims.npg),dtype=np.double,order='c') + double [:] mean = np.empty((Gr.dims.nlg[2],),dtype=np.double,order='c') @@ -240,13 +242,21 @@ cdef class TKE: cdef: Py_ssize_t i Py_ssize_t npg = Gr.dims.npg + double delta = (Gr.dims.dx[0] * Gr.dims.dx[1] * Gr.dims.dx[2])**(1.0/3.0) + double [:] prt = np.zeros((Gr.dims.npg),dtype=np.double,order='c') + double [:] mixing_length = np.zeros((Gr.dims.npg),dtype=np.double,order='c') with nogil: for i in xrange(npg): - prt[i] = DV.values[visc_shift + i]/(DV.values[diff_shift + i] + 1.0e-20) + mixing_length[i] = tke_ell(self.cn, PV.values[e_shift+i], DV.values[bf_shift+i], delta) + prt[i] = delta/(delta + 2.0*mixing_length[i]) mean = Pa.HorizontalMean(Gr,&prt[0]) NS.write_profile('tke_prandtl_number',mean[Gr.dims.gw:-Gr.dims.gw],Pa) + mean = Pa.HorizontalMean(Gr,&mixing_length[0]) + NS.write_profile('tke_mixing_length',mean[Gr.dims.gw:-Gr.dims.gw],Pa) + + return diff --git a/Simulation3d.pyx b/Simulation3d.pyx index 58e9b1cc..9ee46807 100644 --- a/Simulation3d.pyx +++ b/Simulation3d.pyx @@ -158,7 +158,7 @@ class Simulation3d: self.Th.stats_io(self.Gr, self.PV, self.StatsIO, self.Pa) self.Sur.stats_io(self.Gr, self.StatsIO, self.Pa) self.SGS.stats_io(self.Gr,self.DV,self.PV,self.Ke,self.StatsIO,self.Pa) - self.Ke.stats_io(self.Gr,self.PV,self.StatsIO,self.Pa) + self.Ke.stats_io(self.Gr,self.Ref,self.PV,self.StatsIO,self.Pa) return def force_io(self): @@ -171,6 +171,6 @@ class Simulation3d: self.Th.stats_io(self.Gr, self.PV, self.StatsIO, self.Pa) self.Sur.stats_io(self.Gr, self.StatsIO, self.Pa) self.SGS.stats_io(self.Gr,self.DV,self.PV,self.Ke,self.StatsIO,self.Pa) - self.Ke.stats_io(self.Gr,self.PV,self.StatsIO,self.Pa) + self.Ke.stats_io(self.Gr,self.Ref,self.PV,self.StatsIO,self.Pa) return