Alpine3D  Alpine3D-3.1.0
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
SnowDrift.h
Go to the documentation of this file.
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 
41 typedef mio::Array2D<int> CElementArray;
42 typedef mio::Array1D<double> CDoubleArray;
43 typedef mio::Array1D<int> CIntArray;
44 
48 
49 #pragma GCC diagnostic push
50 #pragma GCC diagnostic ignored "-Weffc++"
51 
52 class SnowpackInterface;
53 class EnergyBalance;
54 
55 typedef enum DRIFT_OUTPUT_ENUM {OUT_CONC, OUT_SUBL} DRIFT_OUTPUT;
56 typedef enum PARAM_TYPES {CON,HUM,SUB,TEM,SUB2} param_type;
57 typedef enum ASPECT_TYPES {OTHER,BOTTOM} aspect_type;
58 
59 struct WIND_FIELD {unsigned int start_step;std::string wind;};
60 
71 class SnowDriftA3D {
72  public:
73  SnowDriftA3D(const mio::DEMObject& dem, const mio::Config& cfg);
74 
75  virtual ~SnowDriftA3D();
76 
77  virtual void setSnowSurfaceData(const mio::Grid2DObject& cH_in, const mio::Grid2DObject& sp_in, const mio::Grid2DObject& rg_in,
78  const mio::Grid2DObject& N3_in, const mio::Grid2DObject& rb_in);
79 
80  void setSnowPack(SnowpackInterface &mysnowpack);
81  void setEnergyBalance(EnergyBalance &myeb);
82 
83  virtual void Compute(const mio::Date& calcDate);
84  bool isNewWindField(const unsigned int current_step) /*const*/;
85  void 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,
86  const mio::Grid2DObject& new_rh, const mio::Grid2DObject& new_ta, const mio::Grid2DObject& new_ilwr, const mio::Date& calcDate,
87  const std::vector<mio::MeteoData>& vecMeteo);
88 
89  void GetTResults(double outtime_v[15], double outtime_tau[15], double outtime_salt[15], double outtime_diff[15]);
90 
91  double getTiming() const;
92 
93  void Destroy();
94 
95  std::string getGridsRequirements() const;
96 
97  protected:
98  void Initialize();
99  void ConstructElements();
100  void InitializeNodes(const mio::Grid3DObject& z_readMatr);
101  void CompleteNodes();
102 
103  virtual void compSaltation(bool setbound);
104  virtual void SnowMassChange(bool setbound, const mio::Date& calcDate);
105  virtual void Suspension();
106  virtual void Diffusion(double deltaT, double &diff_max, double t);
107 
108 
109  //sublimation functions begin
110  virtual void Sublimation();
111  double terminalFallVelocity(const double Radius, const double Temperature, const double RH, const double altitude);
112  double calcS(const double concentration, const double sublradius, const double dmdt);
113  double calcSubldM(const double Radius, const double AirTemperature, const double RH,const double WindSpeed, const double altitude);
114  double reynoldsNumberforFallingParticles(const double Radius, const double Windspeed, const double AirTemperature, const double RH, const double altitude);
115  double ventilationVelocity(const double Radius, const double Windspeed,const double AirTemperature, const double RH, const double altitude);
116  double waterVaporDensity(const double Temperature,const double VaporPressure);
117  double RH_from_q(const double AirTemp, const double q, const double altitude);
118  void initializeTRH();
119  //sublimation functions end
120 
121  //---------------------------------------------------------------------
122  //functions which are required for the fe numerics
123  //defined in SnowDriftFENumerics.cc
124  //---------------------------------------------------------------------
125  //bare numerics for element computations, i.e. integrations, jacobian and auxiliary stuff
126  virtual void Jacobian(double *DETERMINANTJ,double J[][3],const int element, const double *P, int k, const int ix, const int iy, const int iz);
127  virtual void J0fun(double J0[3][3], const double J[3][3]);
128  virtual double GQIntB(double *DETERMINANTJ,const int i, const int j);
129  virtual 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]);
130  virtual double GQIntApdx(double DETERMINANTJ[],const double J0M[3][3][8],const int i, const int j, double b[],const double deltak);
131  virtual double GQIntAdxdx(double *DETERMINANTJ, double J0M[3][3][8],const int i, const int j, double *b,const double deltak);
132  virtual void TTfun(double TT[3][8],const double P[]); // P is a pointer pointing at an adress containing a double
133  virtual void phi(double*PHI, double*P);
134  virtual void setQuadraturePoints();
135 
136  //functions required for solving the linear system
137  virtual void matmult(CDoubleArray& res, const CDoubleArray& x, double* sm, int* ijm);
138  virtual void matmult(CDoubleArray& res, const CDoubleArray& x, const CDoubleArray& sA, const CIntArray& colA, CIntArray& rowA);
139  virtual void transmult(CDoubleArray& res, const CDoubleArray& x,double* sm, int* ijm);
140  virtual void SolveEquation(int timeStep, int maxTimeStep, const param_type param );
141  virtual void bicgStab(CDoubleArray& result, CDoubleArray& rhs, const CDoubleArray& sA, const CIntArray& colA, CIntArray& rowA, const int nmax, const double tol, double& testres);
142 
143 
144  //---------------------------------------------------------------------
145  // finite Element functions, basically wrappers for the bare
146  // numerics. Basically two different "classes" of functions, the
147  // first (assembleSystem, applyBoundaryValues, computeDepositionFlux)
148  // contain loops over all or a particular subset of elements and then
149  // invoke the other class of functions which wrap the numerics for
150  // each element (computeDiffusionTensor, computeDriftVector,
151  // computeElementParameter,
152  // computeElementSystem,computeDirichletBoundaryValues,
153  // addElementMatrix
154  //
155  // all these functions are defined in SnowDriftFEControl.cc
156  //---------------------------------------------------------------------
157  virtual void assembleSystem( CIntArray& colA, CIntArray& rowA, CDoubleArray& sA, CDoubleArray& sB, CDoubleArray& Psi, CDoubleArray& f, const double dt);
158  virtual void applyBoundaryValues(CDoubleArray& c00, CDoubleArray& Psi);
159  virtual void prepareSolve();
160 
161 
162  //functions which are required for building the element matrices,
163  virtual void computeDiffusionTensor(double K[3][3], const unsigned int ix, const unsigned int iy, const unsigned int iz);
164  virtual void computeDriftVector(double b[3], const unsigned int ix, const unsigned int iy, const unsigned int iz );
165  virtual 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);
166 
167  virtual 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);
168 
169  virtual 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);
170 
171  virtual 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);
172 
173  virtual void computeDepositionFlux(const CDoubleArray& c, const double theta);
174  virtual void computeDepositionFluxSublimation(const CDoubleArray& c, const double theta);
175 
176 
177  //---------------------------------------------------------------------
178  // Functions for initializing the finite element procedure
179  // defined in SnowDriftFEInit.cc
180  //---------------------------------------------------------------------
181  void setBC_BottomLayer(CDoubleArray& var00, const param_type param);
182  void values_nodes_to_elements(const mio::Grid3DObject& nodesGrid, CDoubleArray& elementsArray );
183  void values_elements_to_nodes(mio::Grid3DObject& nodesGrid, const CDoubleArray& elementsArray );
184  void setRobinBoundaryCondition(const aspect_type aspect, const double gamma_val, const int ix, const int iy, const int iz, CDoubleArray& var00, const param_type param);
185  virtual void InitializeFEData();
186  virtual void prepareSparseMatrix( CIntArray& colA, CIntArray& rowA, CDoubleArray& adjA);
187  virtual void initializeSystem( CIntArray& colA,CIntArray& rowA,CDoubleArray& sA,CDoubleArray& sB, CDoubleArray& rhs,CDoubleArray& f, CDoubleArray& Psi,CDoubleArray& var,CDoubleArray& var00, const param_type param);
188  virtual void resetArray(CDoubleArray& sA);
189  virtual void resetArray(CIntArray& sA);
190  virtual void classifySubdomain();
191  virtual int numberOfNonzeros();
192  void iterativeSublimationCalculation(int timeStep, int maxTimeStep);
193 
194  protected:
195  Saltation saltation_obj;
196 
197  //Time dependent data output after each computation step (SnowDrift::Compute)
198  double time_v[15], time_tau[15], time_salt[15], time_diff[15];
202 
203  mio::IOManager io;
206  mio::Timer timer;
207 
208  mio::Grid2DObject cH;
209  mio::Grid2DObject sp;
210  mio::Grid2DObject rg;
211  mio::Grid2DObject N3;
212  mio::Grid2DObject rb;
213 
214  //the dimensions of the rectangular domain
215  unsigned int nx, ny, nz;
216 
217  // Variables for the finite element solution
218  unsigned int nDOF; //the total number of degrees of freedom
219  unsigned int nNodes; //the total number of nodes
220  unsigned int nElements; //the total number of elements
221  unsigned int nNZ; //the total number of nonzero elements of the system matrix
222 
223  //create system matrix in compressed sparse row (CSR) storage format
228 
229  //auxiliary vector for incorporating inhomogeneous Dirichlet
230  //boundary conditions
232 
233  //vector of nodal concentrations, solution of the linear system
234  CDoubleArray c; //snow concentration in suspension
235  CDoubleArray q; //specific humidity
236  CDoubleArray T; //potential temperature
237 
238  //vector which contains the right hand side of the linear system
240 
241  //vector of sink/source terms
243 
244  //vector which contains boundary and initial conditions
248 
249  //vector which contains boundary and initial conditions
251  //mio::Grid3DObject newElements_precond;
252  //LH_BC
256  void zeroRow(int node);
257 
258  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
259  //LH_DEBUG: Algorithm test BEGIN
262  //LH_DEBUG: Algorithm test END
263  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
264 
265 
266  double qualla;// paramater to vary the deltak coefficient of the SUPG approach
267 
268  double theta; // here you can vary the heigth above the point used in
269 
270  //the face,bar,corner interior arrays are used for the loops over
271  //all elements and are set in classifySubdomain, see there for details
272  mio::Array2D<int> nz_face;
273  mio::Array2D<int> ny_face;
274  mio::Array2D<int> nx_face;
275 
276  mio::Array2D<int> nz_bar;
277  mio::Array2D<int> ny_bar;
278  mio::Array2D<int> nx_bar;
279 
280  mio::Array2D<int> nz_interior;
281  mio::Array2D<int> ny_interior;
282  mio::Array2D<int> nx_interior;
283 
285 
286  //for the quadrature points of the numerical integration scheme
287  mio::Array2D<double> qPoint;
288 
289  //array for mapping the local node indices onto the global nodes
290  //indices (actually this is the analog of the array elems of the old code)
292 
293  //array for mapping the local node indices onto the global
294  //degree-of-freedom (dof) indices
296 
298 
300 
303 
304  //Meteo 2D data
305  mio::Grid2DObject mns, vw, rh, ta, p, psum, psum_ph/*, iswr, ea*/; // TODO ISWR activate, TODO EA activate
306 
307  //Meteo 1D data
308  double ta_1D;
309 
310  mio::Date skip_date; //time step to skip because there would be no drift
311 
312  //3D
314  mio::Grid3DObject nodes_u, nodes_v, nodes_w, nodes_K;
315  mio::Grid3DObject nodes_e, nodes_c;
318 
319  protected:
320  void buildWindFieldsTable(const std::string& wind_field_string);
321  std::vector<struct WIND_FIELD> wind_fields;
323 
324  void debugOutputs(const mio::Date& calcDate, const std::string& fname, const DRIFT_OUTPUT& filetype);
325  void writeOutput(const std::string& fname); //HACK: this should be done by MeteoIO
326 
327  //sublimation constants
329 
330  // constants originally from Snowpack
331  static const double c_red, grain_size, tau_thresh, z0;
332  static const bool thresh_snow;
333 };
334 
335 #pragma GCC diagnostic pop
336 
337 #endif
338 
339 
340 
341 
342 
mio::Grid2DObject vw
Definition: SnowDrift.h:305
virtual void InitializeFEData()
Initializes system size and all data members which are required for the finite element solution...
Definition: SnowDriftFEInit.cc:32
double ta_1D
Definition: SnowDrift.h:308
mio::Grid3DObject nodes_z
Definition: SnowDrift.h:313
void zeroRow(int node)
Definition: SnowDriftFEControl.cc:1029
double time_v[15]
Definition: SnowDrift.h:198
DRIFT_OUTPUT
Definition: SnowDrift.h:55
mio::Grid2DObject sp
Definition: SnowDrift.h:209
Definition: SnowDrift.h:55
mio::Grid2DObject mns
Definition: SnowDrift.h:305
mio::Grid3DObject nodes_e
Definition: SnowDrift.h:315
CDoubleArray T00
Definition: SnowDrift.h:247
mio::Grid3DObject nodes_Subl_ini
Definition: SnowDrift.h:316
double calcS(const double concentration, const double sublradius, const double dmdt)
Calculate sublimation The mass change of a single particle is now extended to all particles...
Definition: Sublimation.cc:100
double ventilationVelocity(const double Radius, const double Windspeed, const double AirTemperature, const double RH, const double altitude)
Ventilation velocity.
Definition: Sublimation.cc:129
Definition: SnowDrift.h:56
double station_altitude
Definition: SnowDrift.h:201
virtual void prepareSparseMatrix(CIntArray &colA, CIntArray &rowA, CDoubleArray &adjA)
prepare sparse matrix Description : generates the auxiliary arrays (i.e. the arrays row-pointer and c...
Definition: SnowDriftFEInit.cc:248
CElementArray elems
Definition: SnowDrift.h:301
void setSnowPack(SnowpackInterface &mysnowpack)
Definition: SnowDrift.cc:707
Definition: SnowDrift.h:57
virtual void assembleSystem(CIntArray &colA, CIntArray &rowA, CDoubleArray &sA, CDoubleArray &sB, CDoubleArray &Psi, CDoubleArray &f, const double dt)
Assemble System loop over all elements, updates the system matrices A, B and the 'bare' right hand si...
Definition: SnowDriftFEControl.cc:43
unsigned int ny
Definition: SnowDrift.h:215
mio::Array2D< double > qPoint
Definition: SnowDrift.h:287
unsigned int nNZ
Definition: SnowDrift.h:221
std::string wind
Definition: SnowDrift.h:59
mio::Array2D< double > c_salt
Definition: SnowDrift.h:302
mio::Grid3DObject nodes_sy
Definition: SnowDrift.h:313
CElementArray nodeMap
Definition: SnowDrift.h:291
CIntArray colA
Definition: SnowDrift.h:226
mio::Grid2DObject p
Definition: SnowDrift.h:305
virtual double GQIntB(double *DETERMINANTJ, const int i, const int j)
GQIntB returns the integral of the product of basis functions i and j over a hexahedral element using...
Definition: SnowDriftFENumerics.cc:170
mio::Grid3DObject nodes_Tair
Definition: SnowDrift.h:316
Definition: SnowDrift.h:59
mio::Grid3DObject nodes_RH
Definition: SnowDrift.h:316
Definition: SnowDrift.h:56
void setBC_BottomLayer(CDoubleArray &var00, const param_type param)
Set Bottom Boundary Condition Set the boundary condition of bottom layer +intialize concentration in ...
Definition: SnowDriftFEInit.cc:320
mio::Grid3DObject nodes_wstar
Definition: SnowDrift.h:313
static const double thermalConductivityofAtm
Definition: SnowDrift.h:328
mio::Array2D< int > nz_face
Definition: SnowDrift.h:272
mio::Grid3DObject nodes_slope
Definition: SnowDrift.h:313
virtual void TTfun(double TT[3][8], const double P[])
TT function TTfun is a function calculating the values of the partial derivative of phi at the point ...
Definition: SnowDriftFENumerics.cc:421
mio::Grid3DObject nodes_K
Definition: SnowDrift.h:314
CDoubleArray c
Definition: SnowDrift.h:234
static const double tau_thresh
Definition: SnowDrift.h:331
virtual void initializeSystem(CIntArray &colA, CIntArray &rowA, CDoubleArray &sA, CDoubleArray &sB, CDoubleArray &rhs, CDoubleArray &f, CDoubleArray &Psi, CDoubleArray &var, CDoubleArray &var00, const param_type param)
Initialize system initializes vectors for the algorithm required for each time step.
Definition: SnowDriftFEInit.cc:390
virtual void setQuadraturePoints()
Definition: SnowDriftFENumerics.cc:50
CDoubleArray gamma
Definition: SnowDrift.h:255
CDoubleArray sB
Definition: SnowDrift.h:225
mio::Grid3DObject nodes_u
Definition: SnowDrift.h:314
virtual ~SnowDriftA3D()
Definition: SnowDrift.cc:160
mio::Grid3DObject nodes_v
Definition: SnowDrift.h:314
mio::Grid2DObject cH
Definition: SnowDrift.h:208
Definition: SnowDrift.h:56
param_type
Definition: SnowDrift.h:56
void InitializeNodes(const mio::Grid3DObject &z_readMatr)
Initialize nodes Each (3D) node receives its (x,y,z), slope, aspect, normals.
Definition: SnowDrift.cc:231
Saltation saltation_obj
Definition: SnowDrift.h:195
CDoubleArray flux_z_subl
Definition: SnowDrift.h:297
double RH_from_q(const double AirTemp, const double q, const double altitude)
Relative humidity Gives the relative humidity for a given air temperature and specific humidity...
Definition: Sublimation.cc:166
virtual void classifySubdomain()
Description... this functions sets summation limits for the loops over boundary elements (faces...
Definition: SnowDriftFEInit.cc:118
void ConstructElements()
Definition: SnowDrift.cc:312
double theta
Definition: SnowDrift.h:268
mio::Array2D< int > nx_face
Definition: SnowDrift.h:274
mio::Grid2DObject rg
Definition: SnowDrift.h:210
virtual void setSnowSurfaceData(const mio::Grid2DObject &cH_in, const mio::Grid2DObject &sp_in, const mio::Grid2DObject &rg_in, const mio::Grid2DObject &N3_in, const mio::Grid2DObject &rb_in)
Definition: SnowDrift.cc:717
virtual void Compute(const mio::Date &calcDate)
Main: Calls the essential routines.
Definition: SnowDrift.cc:655
mio::Array2D< int > nz_bar
Definition: SnowDrift.h:276
virtual void phi(double *PHI, double *P)
Phi - shape functions Calculates the value of the 8 shape functions at point P.
Definition: SnowDriftFENumerics.cc:462
Definition: SnowDrift.h:71
CIntArray n_corner
Definition: SnowDrift.h:284
static const bool thresh_snow
Definition: SnowDrift.h:332
double time_salt[15]
Definition: SnowDrift.h:198
void CompleteNodes()
Initialize concentration according to the rain rate (uniform) Now calculate the missing parameters Sl...
Definition: SnowDrift.cc:170
static const double molecularWeightofWater
Definition: SnowDrift.h:328
Definition: EnergyBalance.h:69
mio::Grid3DObject nodes_w
Definition: SnowDrift.h:314
SnowDriftA3D(const mio::DEMObject &dem, const mio::Config &cfg)
Definition: SnowDrift.cc:50
mio::Grid3DObject nodes_q_ini
Definition: SnowDrift.h:316
std::vector< struct WIND_FIELD > wind_fields
Definition: SnowDrift.h:321
virtual 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)
Compute parameters of the SUPG method Compute the parameters of the SUPG method for each element...
Definition: SnowDriftFEControl.cc:453
int DOITERATION
Definition: SnowDrift.h:199
double auxLayerHeight
Definition: SnowDrift.h:200
Definition: SnowDrift.h:56
double waterVaporDensity(const double Temperature, const double VaporPressure)
CDoubleArray flux_x_subl
Definition: SnowDrift.h:297
Definition: SnowpackInterface.h:121
mio::Grid2DObject ta
Definition: SnowDrift.h:305
double getTiming() const
Definition: SnowDrift.cc:155
virtual void Suspension()
Definition: Suspension.cc:40
virtual void computeDepositionFlux(const CDoubleArray &c, const double theta)
compute deposition flux Authors : Marc Ryser Description : calculates the deposition flux i...
Definition: SnowDriftFEControl.cc:839
CIntArray nnzA
Definition: SnowDrift.h:260
CDoubleArray gNeumann
Definition: SnowDrift.h:253
double reynoldsNumberforFallingParticles(const double Radius, const double Windspeed, const double AirTemperature, const double RH, const double altitude)
Reynolds number Calculate the Reynolds number for a falling particle.
Definition: Sublimation.cc:114
CDoubleArray flux_x
Definition: SnowDrift.h:297
int wind_field_index
Definition: SnowDrift.h:322
static const double grain_size
Definition: SnowDrift.h:331
unsigned int nElements
Definition: SnowDrift.h:220
void initializeTRH()
Initialize RH, T and q Initialize fields of humidity and temperature, necessary for sublimation calcu...
Definition: SnowDrift.cc:838
mio::Array1D< int > CIntArray
Definition: SnowDrift.h:43
virtual void SnowMassChange(bool setbound, const mio::Date &calcDate)
Definition: SnowDrift.cc:390
bool new_wind_status
Definition: SnowDrift.h:299
virtual 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])
GQIntC returns the integral of the product of the gradients of basis functions i and j over a hexahed...
Definition: SnowDriftFENumerics.cc:203
mio::Grid3DObject nodes_Tair_ini
Definition: SnowDrift.h:316
mio::IOManager io
Definition: SnowDrift.h:203
mio::Timer timer
Definition: SnowDrift.h:206
mio::Grid3DObject nodes_WindVel
Definition: SnowDrift.h:316
Definition: SnowDrift.h:56
virtual void applyBoundaryValues(CDoubleArray &c00, CDoubleArray &Psi)
applyBoundaryValues Apply boundary values to all elements Author: Marc Ryser
Definition: SnowDriftFEControl.cc:206
CDoubleArray flux_y_subl
Definition: SnowDrift.h:297
CDoubleArray flux_z
Definition: SnowDrift.h:297
mio::Grid2DObject rh
Definition: SnowDrift.h:305
mio::Grid3DObject nodes_c
Definition: SnowDrift.h:315
CDoubleArray adjA
Definition: SnowDrift.h:261
mio::Array2D< int > nz_interior
Definition: SnowDrift.h:280
CDoubleArray q
Definition: SnowDrift.h:235
virtual void computeDepositionFluxSublimation(const CDoubleArray &c, const double theta)
compute deposition flux Authors : Marc Ryser Description : calculates the deposition flux i...
Definition: SnowDriftFEControl.cc:940
mio::Array2D< int > ny_face
Definition: SnowDrift.h:273
virtual void computeDiffusionTensor(double K[3][3], const unsigned int ix, const unsigned int iy, const unsigned int iz)
Compute diffusion tensor Computes the diffusion tensor of an element, accesses the nodes array and th...
Definition: SnowDriftFEControl.cc:343
void values_elements_to_nodes(mio::Grid3DObject &nodesGrid, const CDoubleArray &elementsArray)
elements_to_nodes Extract a 3D-grid of nodes from an array of elements
Definition: SnowDriftFEInit.cc:511
mio::Array2D< int > ny_interior
Definition: SnowDrift.h:281
mio::Array2D< double > mns_nosubl
Definition: SnowDrift.h:302
unsigned int nx
Definition: SnowDrift.h:215
CDoubleArray precond
Definition: SnowDrift.h:250
void 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, const mio::Grid2DObject &new_rh, const mio::Grid2DObject &new_ta, const mio::Grid2DObject &new_ilwr, const mio::Date &calcDate, const std::vector< mio::MeteoData > &vecMeteo)
Sets the required meteo fields.
Definition: SnowDrift.cc:780
virtual void Sublimation()
Sublimation Calculates drifting snow sublimation at every node of the concentration field...
Definition: Sublimation.cc:33
void Initialize()
Initialize for snowdrift Allocate saltation variable, snow height and new snow mass per bottom elemen...
Definition: SnowDrift.cc:122
mio::Grid3DObject nodes_Subl
Definition: SnowDrift.h:316
mio::Grid3DObject nodes_q
Definition: SnowDrift.h:316
virtual void transmult(CDoubleArray &res, const CDoubleArray &x, double *sm, int *ijm)
transmult this is the transmult function that multiplies the vector x by the transpose of the matrix ...
Definition: SnowDriftFENumerics.cc:536
virtual void Diffusion(double deltaT, double &diff_max, double t)
Definition: Suspension.cc:31
mio::Array2D< double > saltation
Definition: SnowDrift.h:302
virtual void compSaltation(bool setbound)
Calculate Saltation Fluxes for all Bottom Elements.
Definition: Saltation.cc:51
mio::Array2D< int > nx_bar
Definition: SnowDrift.h:278
mio::Array1D< double > CDoubleArray
Definition: SnowDrift.h:42
virtual void computeDriftVector(double b[3], const unsigned int ix, const unsigned int iy, const unsigned int iz)
Compute the drift vector of an element.
Definition: SnowDriftFEControl.cc:421
void Destroy()
Definition: SnowDrift.cc:153
virtual double GQIntAdxdx(double *DETERMINANTJ, double J0M[3][3][8], const int i, const int j, double *b, const double deltak)
GQIntAdxdx ...
Definition: SnowDriftFENumerics.cc:351
static const double c_red
Definition: SnowDrift.h:331
double qualla
Definition: SnowDrift.h:266
virtual double GQIntApdx(double DETERMINANTJ[], const double J0M[3][3][8], const int i, const int j, double b[], const double deltak)
GQIntApdx returns the integral of the product of a basis function i and the gradient of basis functio...
Definition: SnowDriftFENumerics.cc:297
mio::Grid2DObject psum
Definition: SnowDrift.h:305
Definition: SnowDrift.h:55
unsigned int nNodes
Definition: SnowDrift.h:219
Definition: SnowDrift.h:57
virtual 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)
compute element system ....
Definition: SnowDriftFEControl.cc:540
virtual int numberOfNonzeros()
Number of nonzero elements of the finite element system matrix Description : computes the number of n...
Definition: SnowDriftFEInit.cc:216
double time_diff[15]
Definition: SnowDrift.h:198
double calcSubldM(const double Radius, const double AirTemperature, const double RH, const double WindSpeed, const double altitude)
Calculate sublimation mass change This function calculates the mass change of a single ice particle f...
Definition: Sublimation.cc:60
void setEnergyBalance(EnergyBalance &myeb)
Definition: SnowDrift.cc:712
aspect_type
Definition: SnowDrift.h:57
void buildWindFieldsTable(const std::string &wind_field_string)
Definition: SnowDrift.cc:81
static int steps
Definition: AlpineMain.cc:47
CDoubleArray Psi
Definition: SnowDrift.h:231
CElementArray dofMap
Definition: SnowDrift.h:295
bool STATIONARY
Definition: SnowDrift.h:317
CDoubleArray T
Definition: SnowDrift.h:236
void setRobinBoundaryCondition(const aspect_type aspect, const double gamma_val, const int ix, const int iy, const int iz, CDoubleArray &var00, const param_type param)
Set Robin Boundary Condition.
Definition: SnowDriftFEInit.cc:362
virtual void SolveEquation(int timeStep, int maxTimeStep, const param_type param)
Solve the advection diffusion equation.
Definition: Suspension.cc:123
SnowpackInterface * snowpack
Definition: SnowDrift.h:204
static const double z0
Definition: SnowDrift.h:331
CDoubleArray c00
Definition: SnowDrift.h:245
virtual void J0fun(double J0[3][3], const double J[3][3])
J0 function computes the matrix of cofactors of the input matrix J and stores it in J0...
Definition: SnowDriftFENumerics.cc:144
CDoubleArray sA
Definition: SnowDrift.h:224
unsigned int nz
Definition: SnowDrift.h:215
std::string getGridsRequirements() const
Definition: SnowDrift.cc:76
CDoubleArray gDirichlet
Definition: SnowDrift.h:254
virtual void bicgStab(CDoubleArray &result, CDoubleArray &rhs, const CDoubleArray &sA, const CIntArray &colA, CIntArray &rowA, const int nmax, const double tol, double &testres)
bicgStab iterative equation solver iterative equation solver Tests : Tested by the followin procedure...
Definition: SnowDriftFENumerics.cc:569
CDoubleArray rhs
Definition: SnowDrift.h:239
bool isNewWindField(const unsigned int current_step)
Definition: SnowDrift.cc:104
CDoubleArray f
Definition: SnowDrift.h:242
double time_tau[15]
Definition: SnowDrift.h:198
unsigned int start_step
Definition: SnowDrift.h:59
void iterativeSublimationCalculation(int timeStep, int maxTimeStep)
Calculate the steady state sublimation in several steps. Single feedbacks can be switched on or off i...
Definition: Suspension.cc:199
mio::Array2D< int > CElementArray
Definition: SnowDrift.h:41
mio::Array2D< double > dif_mns_subl
Definition: SnowDrift.h:302
mio::Array2D< double > mns_subl
Definition: SnowDrift.h:302
EnergyBalance * eb
Definition: SnowDrift.h:205
CDoubleArray q00
Definition: SnowDrift.h:246
mio::Array2D< int > ny_bar
Definition: SnowDrift.h:277
virtual void prepareSolve()
prepareSolve Some preparations to solve a 3D-field.
Definition: Suspension.cc:361
virtual 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)
Add element matrix given a sparse matrix A in CSR format ( specified by sA, colInd, rowPtr) and the element matrix Bel, this function adds the element contributions to the global ones ; the first length_spec entries of the vector spec contains the indices of the matrix which are inserted Comments : Beware of the enumeration of the matrix sA (size is one bigger than required, zero index is unused => dangerous accesses nodeMap, nodeMap.
Definition: SnowDriftFEControl.cc:642
mio::Array2D< int > nx_interior
Definition: SnowDrift.h:282
virtual void Jacobian(double *DETERMINANTJ, double J[][3], const int element, const double *P, int k, const int ix, const int iy, const int iz)
Compute Jacobian This function calculates the Jacobian at the given point P and and stores the value ...
Definition: SnowDriftFENumerics.cc:78
double terminalFallVelocity(const double Radius, const double Temperature, const double RH, const double altitude)
Terminal fall velocity NB assumed to be 0.5.
Definition: Sublimation.cc:146
mio::Grid2DObject N3
Definition: SnowDrift.h:211
void debugOutputs(const mio::Date &calcDate, const std::string &fname, const DRIFT_OUTPUT &filetype)
Definition: SnowDrift.cc:727
static const double USTAR
Definition: SnowDrift.h:328
CIntArray rowA
Definition: SnowDrift.h:227
void writeOutput(const std::string &fname)
Write output Writes the values of several nodes-fields.
Definition: SnowDrift.cc:741
mio::Grid2DObject rb
Definition: SnowDrift.h:212
static const double kinematicViscosityAir
Definition: SnowDrift.h:328
unsigned int nDOF
Definition: SnowDrift.h:218
void GetTResults(double outtime_v[15], double outtime_tau[15], double outtime_salt[15], double outtime_diff[15])
Definition: SnowDrift.cc:698
mio::Date skip_date
Definition: SnowDrift.h:310
mio::Grid2DObject psum_ph
Definition: SnowDrift.h:305
mio::Grid3DObject nodes_tmp_c
Definition: SnowDrift.h:316
mio::Grid3DObject nodes_y
Definition: SnowDrift.h:313
void values_nodes_to_elements(const mio::Grid3DObject &nodesGrid, CDoubleArray &elementsArray)
nodes_to_elements
Definition: SnowDriftFEInit.cc:473
mio::Grid3DObject nodes_x
Definition: SnowDrift.h:313
virtual 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)
Definition: SnowDriftFEControl.cc:721
virtual void matmult(CDoubleArray &res, const CDoubleArray &x, double *sm, int *ijm)
matmult computes a matrix vector product for a sparse matrix
Definition: SnowDriftFENumerics.cc:484
virtual void resetArray(CDoubleArray &sA)
Definition: SnowDriftFEInit.cc:461
mio::Grid3DObject nodes_sx
Definition: SnowDrift.h:313
CDoubleArray flux_y
Definition: SnowDrift.h:297