Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/snowdrift/SnowDrift.h

1/***********************************************************************************/
2/* Copyright 2009-2015 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
3/***********************************************************************************/
4/* This file is part of Alpine3D.
5 Alpine3D is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 Alpine3D is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with Alpine3D. If not, see <http://www.gnu.org/licenses/>.
17*/
18#ifndef SNOWDRIFTA3D_H
19#define SNOWDRIFTA3D_H
20
21#define WS0 0.5
22#define TKE 0
23#define SALTATION 1// switch for saltation simulation
24#define SUBLIMATION 0 // switch for drifting snow sublimation
25#define FIELD3D_OUTPUT 0// output with all three-dimensional fields (only possible for sublimation)
26#define SUBLIMATION_OUTPUT 0// debug output of drifting snow sublimation
27#define T_FB 1 //switch for feedback between sublimation and air temperature
28#define Q_FB 1 //switch for feedback between sublimation and humidity
29#define C_FB 1 //switch for feedback between sublimation and snow concentration
30#define READK 0 //define as 1 if you have K from ARPS wind fields INCLUDING turbulence
31#define WRITE_DRIFT_FLUXES 0 //set to 1 in order to write snow drift fluxes
32#define dt_diff 0.5 /* Small calculation step length for snow diffusion */
33
34#include <stdio.h>
35#include <math.h>
36#include <stdlib.h>
37#include <iostream>
38#include <meteoio/MeteoIO.h>
39#include <meteoio/plugins/ARPSIO.h>
40
41typedef mio::Array2D<int> CElementArray;
42typedef mio::Array1D<double> CDoubleArray;
43typedef mio::Array1D<int> CIntArray;
44
45#include <alpine3d/SnowpackInterface.h>
46#include <alpine3d/ebalance/EnergyBalance.h>
47#include <alpine3d/snowdrift/checksum.h>
48
49#pragma GCC diagnostic push
50#pragma GCC diagnostic ignored "-Weffc++"
51
52class SnowpackInterface;
53class EnergyBalance;
54
55typedef enum DRIFT_OUTPUT_ENUM {OUT_CONC, OUT_SUBL} DRIFT_OUTPUT;
56typedef enum PARAM_TYPES {CON,HUM,SUB,TEM,SUB2} param_type;
57typedef enum ASPECT_TYPES {OTHER,BOTTOM} aspect_type;
58
59struct WIND_FIELD {unsigned int start_step;std::string wind;};
60
61/**
62 * @page snowdrift Snowdrift
63 * This module computes the preferential deposition and redistribution of snow by the wind (see \ref principles_snowdrift). It is
64 * enabled using the "--enable-drift" command line option.
65 * The 3D wind fields are read using a GRID3D <A HREF="https://models.slf.ch/p/meteoio">MeteoIO</A> plugin such as ARPS:
66 * @code
67 * GRID3D= ARPS
68 * GRID3DPATH= ../input/wind_fields/
69 * @endcode
70 *
71 * The WINDFIELDS key must be defined in the [Input] section and contains a space delimited list of wind fields files (here within GRID3DPATH) and
72 * associated number of time steps:
73 * @code
74 * WINDFIELDS = sw3.asc 1 nw3.asc 3 ww0.asc 2 nw9.asc 5 nw6.asc 10 ww0.asc 5 sw3.asc 6 nw3.asc 1
75 * @endcode
76 *
77 */
78class SnowDriftA3D {
79public:
80SnowDriftA3D(const mio::DEMObject& dem, const mio::Config& cfg);
81
82virtual ~SnowDriftA3D();
83
84virtual void setSnowSurfaceData(const mio::Grid2DObject& cH_in, const mio::Grid2DObject& sp_in, const mio::Grid2DObject& rg_in,
85 const mio::Grid2DObject& N3_in, const mio::Grid2DObject& rb_in);
86
87void setSnowPack(SnowpackInterface &mysnowpack);
88void setEnergyBalance(EnergyBalance &myeb);
89
90virtual void Compute(const mio::Date& calcDate);
91bool isNewWindField(const unsigned int current_step) /*const*/;
92void setMeteo (const unsigned int& steps, const mio::Grid2DObject& new_psum, const mio::Grid2DObject& new_psum_ph, const mio::Grid2DObject& new_p, const mio::Grid2DObject& new_vw,
93 const mio::Grid2DObject& new_rh, const mio::Grid2DObject& new_ta, const mio::Grid2DObject& new_ilwr, const mio::Date& calcDate,
94 const std::vector<mio::MeteoData>& vecMeteo);
95
96void GetTResults(double outtime_v[15], double outtime_tau[15], double outtime_salt[15], double outtime_diff[15]);
97
98double getTiming() const;
99
100void Destroy();
101
102std::string getGridsRequirements() const;
103
104protected:
105void Initialize();
106void ConstructElements();
107void InitializeNodes(const mio::Grid3DObject& z_readMatr);
108void CompleteNodes();
109
110virtual void compSaltation(bool setbound);
111virtual void SnowMassChange(bool setbound, const mio::Date& calcDate);
112virtual void Suspension();
113virtual void Diffusion(double deltaT, double &diff_max, double t);
114
115
116//sublimation functions begin
117virtual void Sublimation();
118double terminalFallVelocity(const double Radius, const double Temperature, const double RH, const double altitude);
119double calcS(const double concentration, const double sublradius, const double dmdt);
120double calcSubldM(const double Radius, const double AirTemperature, const double RH,const double WindSpeed, const double altitude);
121double reynoldsNumberforFallingParticles(const double Radius, const double Windspeed, const double AirTemperature, const double RH, const double altitude);
122double ventilationVelocity(const double Radius, const double Windspeed,const double AirTemperature, const double RH, const double altitude);
123double waterVaporDensity(const double Temperature,const double VaporPressure);
124double RH_from_q(const double AirTemp, const double q, const double altitude);
125void initializeTRH();
126//sublimation functions end
127
128//---------------------------------------------------------------------
129//functions which are required for the fe numerics
130//defined in SnowDriftFENumerics.cc
131//---------------------------------------------------------------------
132//bare numerics for element computations, i.e. integrations, jacobian and auxiliary stuff
133virtual void Jacobian(double *DETERMINANTJ,double J[][3],const int element, const double *P, int k, const int ix, const int iy, const int iz);
134virtual void J0fun(double J0[3][3], const double J[3][3]);
135virtual double GQIntB(double *DETERMINANTJ,const int i, const int j);
136virtual double GQIntC(double * DETERMINANTJ, const double J0M[3][3][8], const int i, const int j, const double b[3],const double K[3][3]);
137virtual double GQIntApdx(double DETERMINANTJ[],const double J0M[3][3][8],const int i, const int j, double b[],const double deltak);
138virtual double GQIntAdxdx(double *DETERMINANTJ, double J0M[3][3][8],const int i, const int j, double *b,const double deltak);
139virtual void TTfun(double TT[3][8],const double P[]);// P is a pointer pointing at an adress containing a double
140virtual void phi(double*PHI, double*P);
141virtual void setQuadraturePoints();
142
143//functions required for solving the linear system
144virtual void matmult(CDoubleArray& res, const CDoubleArray& x, double* sm, int* ijm);
145virtual void matmult(CDoubleArray& res, const CDoubleArray& x, const CDoubleArray& sA, const CIntArray& colA, CIntArray& rowA);
146virtual void transmult(CDoubleArray& res, const CDoubleArray& x,double* sm, int* ijm);
147virtual void SolveEquation(int timeStep, int maxTimeStep, const param_type param );
148virtual void bicgStab(CDoubleArray& result, CDoubleArray& rhs, const CDoubleArray& sA, const CIntArray& colA, CIntArray& rowA, const int nmax, const double tol, double& testres);
149
150
151//---------------------------------------------------------------------
152// finite Element functions, basically wrappers for the bare
153// numerics. Basically two different "classes" of functions, the
154// first (assembleSystem, applyBoundaryValues, computeDepositionFlux)
155// contain loops over all or a particular subset of elements and then
156// invoke the other class of functions which wrap the numerics for
157// each element (computeDiffusionTensor, computeDriftVector,
158// computeElementParameter,
159// computeElementSystem,computeDirichletBoundaryValues,
160// addElementMatrix
161//
162// all these functions are defined in SnowDriftFEControl.cc
163//---------------------------------------------------------------------
164virtual void assembleSystem( CIntArray& colA, CIntArray& rowA, CDoubleArray& sA, CDoubleArray& sB, CDoubleArray& Psi, CDoubleArray& f, const double dt);
165virtual void applyBoundaryValues(CDoubleArray& c00, CDoubleArray& Psi);
166virtual void prepareSolve();
167
168
169//functions which are required for building the element matrices,
170virtual void computeDiffusionTensor(double K[3][3], const unsigned int ix, const unsigned int iy, const unsigned int iz);
171virtual void computeDriftVector(double b[3], const unsigned int ix, const unsigned int iy, const unsigned int iz );
172virtual void computeElementParameters(const int& element, double DETERMINANTJ[8], double J0M[3][3][8], double J0[3][3], double J[3][3], double b[3], double K[3][3], double& deltak, double& qualla, const int ix, const int iy, const int iz);
173
174virtual void computeElementSystem(int &element, int &nDofNodes, int* dofNode, double Ael[9][9], double Del[9][9], bool stationary, double DETERMINANTJ[8], double J0M[3][3][8], double b[3], double K[3][3], double &deltak, const double &dt, CDoubleArray& f, CDoubleArray& Psi);
175
176virtual void computeDirichletBoundaryValues(int element,double DETERMINANTJ[8],double J0M[3][3][8], double J0[3][3], double b[3], double K[3][3], double deltak, int spec[8], int length_spec, int length_complSpec, CDoubleArray& c00, CDoubleArray& Psi);
177
178virtual void addElementMatrix( CDoubleArray& sA, const CIntArray& colInd, const CIntArray& rowPtr,const double Bel[9][9], const int element, const int* spec,const int length_spec);
179
180virtual void computeDepositionFlux(const CDoubleArray& c, const double theta);
181virtual void computeDepositionFluxSublimation(const CDoubleArray& c, const double theta);
182
183
184//---------------------------------------------------------------------
185// Functions for initializing the finite element procedure
186// defined in SnowDriftFEInit.cc
187//---------------------------------------------------------------------
188void setBC_BottomLayer(CDoubleArray& var00, const param_type param);
189void values_nodes_to_elements(const mio::Grid3DObject& nodesGrid, CDoubleArray& elementsArray );
190void values_elements_to_nodes(mio::Grid3DObject& nodesGrid, const CDoubleArray& elementsArray );
191void setRobinBoundaryCondition(const aspect_type aspect, const double gamma_val, const int ix, const int iy, const int iz, CDoubleArray& var00, const param_type param);
192virtual void InitializeFEData();
193virtual void prepareSparseMatrix( CIntArray& colA, CIntArray& rowA, CDoubleArray& adjA);
194virtual void initializeSystem( CIntArray& colA,CIntArray& rowA,CDoubleArray& sA,CDoubleArray& sB, CDoubleArray& rhs,CDoubleArray& f, CDoubleArray& Psi,CDoubleArray& var,CDoubleArray& var00, const param_type param);
195virtual void resetArray(CDoubleArray& sA);
196virtual void resetArray(CIntArray& sA);
197virtual void classifySubdomain();
198virtual int numberOfNonzeros();
199void iterativeSublimationCalculation(int timeStep, int maxTimeStep);
200
201protected:
202Saltation saltation_obj;
203
204//Time dependent data output after each computation step (SnowDrift::Compute)
205double time_v[15], time_tau[15], time_salt[15], time_diff[15];
206int DOITERATION;
207double auxLayerHeight;
208double station_altitude;
209
210mio::IOManager io;
211SnowpackInterface *snowpack;
212EnergyBalance *eb;
213mio::Timer timer;
214
215mio::Grid2DObject cH;
216mio::Grid2DObject sp;
217mio::Grid2DObject rg;
218mio::Grid2DObject N3;
219mio::Grid2DObject rb;
220
221//the dimensions of the rectangular domain
222unsigned int nx, ny, nz;
223
224// Variables for the finite element solution
225unsigned int nDOF; //the total number of degrees of freedom
226unsigned int nNodes; //the total number of nodes
227unsigned int nElements; //the total number of elements
228unsigned int nNZ; //the total number of nonzero elements of the system matrix
229
230//create system matrix in compressed sparse row (CSR) storage format
231CDoubleArray sA;
232CDoubleArray sB;
233CIntArray colA;
234CIntArray rowA;
235
236//auxiliary vector for incorporating inhomogeneous Dirichlet
237//boundary conditions
238CDoubleArray Psi;
239
240//vector of nodal concentrations, solution of the linear system
241CDoubleArray c; //snow concentration in suspension
242CDoubleArray q; //specific humidity
243CDoubleArray T; //potential temperature
244
245//vector which contains the right hand side of the linear system
246CDoubleArray rhs;
247
248//vector of sink/source terms
249CDoubleArray f;
250
251//vector which contains boundary and initial conditions
252CDoubleArray c00;
253CDoubleArray q00;
254CDoubleArray T00;
255
256//vector which contains boundary and initial conditions
257CDoubleArray precond;
258//mio::Grid3DObject newElements_precond;
259//LH_BC
260CDoubleArray gNeumann;
261CDoubleArray gDirichlet;
262CDoubleArray gamma;
263void zeroRow(int node);
264
265//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
266//LH_DEBUG: Algorithm test BEGIN
267CIntArray nnzA;
268CDoubleArray adjA;
269//LH_DEBUG: Algorithm test END
270//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
271
272
273double qualla;// paramater to vary the deltak coefficient of the SUPG approach
274
275double theta; // here you can vary the heigth above the point used in
276
277//the face,bar,corner interior arrays are used for the loops over
278//all elements and are set in classifySubdomain, see there for details
279mio::Array2D<int> nz_face;
280mio::Array2D<int> ny_face;
281mio::Array2D<int> nx_face;
282
283mio::Array2D<int> nz_bar;
284mio::Array2D<int> ny_bar;
285mio::Array2D<int> nx_bar;
286
287mio::Array2D<int> nz_interior;
288mio::Array2D<int> ny_interior;
289mio::Array2D<int> nx_interior;
290
291CIntArray n_corner;
292
293//for the quadrature points of the numerical integration scheme
294mio::Array2D<double> qPoint;
295
296//array for mapping the local node indices onto the global nodes
297//indices (actually this is the analog of the array elems of the old code)
298CElementArray nodeMap;
299
300//array for mapping the local node indices onto the global
301//degree-of-freedom (dof) indices
302CElementArray dofMap;
303
304CDoubleArray flux_x, flux_y, flux_z, flux_x_subl, flux_y_subl, flux_z_subl;
305
306bool new_wind_status;
307
308CElementArray elems;
309mio::Array2D<double> saltation, c_salt, mns_subl, mns_nosubl,dif_mns_subl;
310
311//Meteo 2D data
312mio::Grid2DObject mns, vw, rh, ta, p, psum, psum_ph/*, iswr, ea*/; // TODO ISWR activate, TODO EA activate
313
314//Meteo 1D data
315double ta_1D;
316
317mio::Date skip_date; //time step to skip because there would be no drift
318
319//3D
320mio::Grid3DObject nodes_x, nodes_y, nodes_z, nodes_sy, nodes_sx, nodes_slope, nodes_wstar;
321mio::Grid3DObject nodes_u, nodes_v, nodes_w, nodes_K;
322mio::Grid3DObject nodes_e, nodes_c;
323mio::Grid3DObject nodes_Tair,nodes_Tair_ini, nodes_q, nodes_q_ini, nodes_RH, nodes_Subl, nodes_Subl_ini, nodes_WindVel, nodes_tmp_c;
324bool STATIONARY;
325
326protected:
327void buildWindFieldsTable(const std::string& wind_field_string);
328std::vector<struct WIND_FIELD> wind_fields;
329int wind_field_index;
330
331void debugOutputs(const mio::Date& calcDate, const std::string& fname, const DRIFT_OUTPUT& filetype);
332void writeOutput(const std::string& fname); //HACK: this should be done by MeteoIO
333
334//sublimation constants
335static const double kinematicViscosityAir, USTAR, molecularWeightofWater, thermalConductivityofAtm;
336
337// constants originally from Snowpack
338static const double c_red, grain_size, tau_thresh, z0;
339static const bool thresh_snow;
340};
341
342#pragma GCC diagnostic pop
343
344#endif
345
346
347
348
349

Archive Download this file

Revision: HEAD