Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/snowdrift/SnowDriftFEControl.cc

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#include <assert.h>
19#include <iostream>
20#include <cstdlib>
21#include <string>
22
23#include <alpine3d/snowdrift/Cell.h>
24#include <meteoio/MeteoIO.h>
25
26#include <alpine3d/snowdrift/SnowDrift.h>
27
28using namespace mio;
29
30/**
31 *@brief Assemble System
32 * loop over all elements, updates the system matrices A, B and the 'bare' right hand side (rhs)
33 * of the system, in other words: prepares the system prior to the inclusion of dirichlet
34 * boundary conditions
35 * @param colA column index to locate within sparse matrix
36 * @param rowA row index
37 * @param sA "system matrix"
38 * @param sB "system matrix"
39 * @param Psi vector for incorporating inhomogeneous Dirichlet BC
40 * @param f source in diffusion equation
41 * @param dt
42*/
43void SnowDriftA3D::assembleSystem( CIntArray& colA_loc,
44CIntArray& rowA_loc,
45CDoubleArray& sA_loc,
46CDoubleArray& sB_loc,
47CDoubleArray& Psi_loc,
48CDoubleArray& f_loc,
49const double dt)
50{
51 //-------------------------------------------------------------
52 //variables for integrals over elements
53 //-------------------------------------------------------------
54 double DETERMINANTJ[8]; // to store the determinant of the
55 // isoparametric point transf. at the 8
56 // int. points
57 double J0M[3][3][8]; // to store the eight J0 matrices at the 8
58 // int. points
59 double J[3][3]; // the Jacobian matrix
60 double J0[3][3]; // the J0 matrix
61 //element variables
62 double b[3];// the wind field
63 double K[3][3];// the diffusion matrix
64 //the element matrix
65 double Ael[9][9];
66 double Del[9][9];
67
68 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69 // ELEMENT MATRICES AND ASSEMBLING
70 // now run over all elements: first interior of domain, then faces, corners
71 // and bars
72 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73
74 double deltak; //the SUPG parameter
75 int element; //used as element counter
76
77 //the array dofNode contains the local node indices the first
78 //nDofNode entries are the degrees of freedom of that element and
79 //the following nBoundaryNodes = 8-nDofNodes entries are the
80 //boundary nodes. the number of boundary and dof-nodes depends on
81 //the location of the element (face,bar,corner) and is set by the Cell class
82 int dofNode[8];
83 int nDofNodes;
84 int nBoundaryNodes;
85
86 Cell cell;
87 std::string cellType;
88
89 //LH_BC
90 //loop over all elements
91for ( unsigned int iz = 0; iz < nz-1;iz++){
92for ( unsigned int iy = 0; iy < ny-1;iy++){
93for ( unsigned int ix = 0; ix < nx-1;ix++){
94
95element = iz*(nx-1)*(ny-1)+iy*(nx-1) + ix;
96 cellType = "interior";
97 cell.classifyNodes(dofNode, &nDofNodes, &nBoundaryNodes, cellType,0);
98
99computeDriftVector(b,ix,iy,iz);
100
101computeDiffusionTensor(K,ix,iy,iz);
102
103computeElementParameters(element, DETERMINANTJ,J0M,J0,J,b,K,deltak,qualla,ix,iy,iz);
104
105computeElementSystem(element,nDofNodes,dofNode,Ael,Del,STATIONARY, DETERMINANTJ,J0M,b,K,deltak,dt,f_loc,Psi_loc);
106 addElementMatrix(sA_loc, colA_loc, rowA_loc, Ael,element,dofNode,nDofNodes);//LH
107for ( int i = 0; i < 8; i++){
108 precond[ nodeMap[element][i] ] += Ael[i+1][i+1];
109}
110 addElementMatrix(sB_loc, colA_loc, rowA_loc, Del,element,dofNode,nDofNodes);//LH
111
112 //boundary conditions
113int boundaryFaceNodes[6][4]=
114{
115{0, 1, 4, 5},//neg y direction in global system
116{2, 3, 6, 7},//pos y direction in global system
117{1, 2, 5, 6},//pos x direction in global system
118{0, 3, 4, 7},//neg x direction in global system
119{4, 5, 6, 7},//pos z direction in global system
120{0, 1, 2, 3} //neg z direction in global system
121};
122int isBFace[6];
123isBFace[0] = (iy == 0) ? 1 : 0;
124isBFace[1] = (iy == ny-2)? -1 : 0;
125isBFace[2] = (ix == nx-2)? 1 : 0;
126isBFace[3] = (ix == 0) ? -1 : 0;
127isBFace[4] = (iz == nz-2)? 1 : 0;
128isBFace[5] = (iz == 0) ? -1 : 0;
129
130double BCel[9][9];
131double surfaceMetric;
132double qp[3];
133double PHI[8];
134double rhsel[8];
135
136for ( int i = 0; i < 8; i++){
137rhsel[i] = 0;
138for ( int j = 0; j < 8; j++){
139BCel[i+1][j+1] = 0;
140}
141}
142
143//loop over element faces
144for ( unsigned int bf = 0; bf < 6; bf++) {
145//if b-face
146if ( isBFace[bf] != 0 ){
147//loop over quadrature points
148//compute face averages
149double gc = 0;
150double gN = 0;
151double gD = 0;
152for ( unsigned int k = 0; k < 4; k++ ) {
153gc += 0.25*gamma[ nodeMap[element][ boundaryFaceNodes[bf][k] ] ];
154gN += 0.25*gNeumann[ nodeMap[element][ boundaryFaceNodes[bf][k] ] ];
155gD += 0.25*gDirichlet[ nodeMap[element][ boundaryFaceNodes[bf][k] ] ];
156}
157 //compute surface integral
158 //loop over quadrature points
159for ( unsigned int k = 0; k < 4; k++ ) {
160qp[0]= qPoint(0, boundaryFaceNodes[bf][k]);
161qp[1]= qPoint(1, boundaryFaceNodes[bf][k]);
162qp[2]= qPoint(2, boundaryFaceNodes[bf][k]);
163
164//coordinate direction
165int cDir = bf/2;
166qp[ cDir ] = isBFace[bf];
167phi(PHI,qp);
168Jacobian(DETERMINANTJ,J,element,qp,k,ix,iy,iz);
169J0fun(J0,J);
170
171surfaceMetric = 0;
172 for ( unsigned int l = 0; l < 3; l++) {
173 surfaceMetric += J0[l][cDir]*J0[l][cDir];
174 }
175 surfaceMetric = sqrt(surfaceMetric);
176
177 for ( unsigned int i = 0; i < 8; i++){
178
179rhsel[i] += surfaceMetric*(gc*gD-gN)*PHI[i];
180 for ( unsigned int j = 0; j < 8; j++){
181 BCel[i+1][j+1] += ( surfaceMetric * gc* PHI[i]*PHI[j]);
182}
183 }
184}
185 } //if b-face
186}//loop element faces
187
188for ( unsigned int i = 0; i < 8; i++){
189rhs[ nodeMap[element][i] ] += rhsel[i];
190precond[ nodeMap[element][i] ] += BCel[i+1][i+1];
191}
192
193addElementMatrix(sA_loc,colA_loc,rowA_loc,BCel,element,dofNode,nDofNodes);//LH
194}//end of ix
195}// end of iy
196}//end of iz
197}
198
199/**
200 *@brief applyBoundaryValues
201 * Apply boundary values to all elements
202 * Author: Marc Ryser
203 * @param var00 initial variable
204 * @param Psi vector for incorporating inhomogeneous Dirichlet BC
205*/
206void SnowDriftA3D::applyBoundaryValues(CDoubleArray& var00,
207 CDoubleArray& Psi_loc)
208{
209 //variables for integrals over elements
210 double DETERMINANTJ[8]; // to store the determinant of the
211 // isoparametric point transf. at the 8
212 // int. points
213 double J0M[3][3][8]; // to store the eight J0 matrices at the 8
214 // int. points
215 double J[3][3]; // the Jacobian matrix
216 double J0[3][3]; // the J0 matrix
217 //element variables
218 double b[3];// the wind field
219 double K[3][3];// the diffusion matrix
220 double deltak; //the SUPG parameter
221 int element; //used as element counter
222
223 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
224 // Loop over all boundary elements
225 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
226 //the array dofNode contains the local node indices the first
227 //nDofNode entries are the degrees of freedom of that element and
228 //the following nBoundaryNodes = 8-nDofNodes entries are the
229 //boundary nodes. the number of boundary and dof-nodes depends on
230 //the location of the element (face,bar,corner) and is set by the Cell class
231 int dofNode[8];
232 int nDofNodes;
233 int nBoundaryNodes;
234
235 Cell cell;
236 std::string cellType;
237
238 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
239 // Step 1: ELEMENTS AND ELEMENT MATRICES ON FACES
240 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
241 size_t nz_face_size, tmp;
242 nz_face.size(nz_face_size, tmp);
243 //for (int face = 0; face < nz_face.size(); face++)
244 for (size_t face = 0; face < nz_face_size; face++) {
245for (int layer = nz_face[face][0]; layer<=nz_face[face][1]; layer++ ){
246 for (int depth = ny_face[face][0]; depth<=ny_face[face][1]; depth++ ){
247for (int width = nx_face[face][0]; width<=nx_face[face][1]; width++ ){
248element = layer*(nx-1)*(ny-1)+depth*(nx-1)+width;
249computeDriftVector(b,width,depth,layer);
250computeDiffusionTensor(K,width,depth,layer);
251cellType = "face";
252cell.classifyNodes(dofNode,&nDofNodes,&nBoundaryNodes,cellType,(int)face);
253computeElementParameters(element, DETERMINANTJ,J0M,J0,J,b,K,deltak,qualla,width,depth,layer);
254
255//DirichletBoundaryValues hasn't been implemented
256computeDirichletBoundaryValues( element,DETERMINANTJ,J0M,J0,b,K,deltak,dofNode,nDofNodes,nBoundaryNodes,var00,Psi_loc);
257} // end of width loop
258 }// end of depth loop
259}// end of layer loop
260 }// end of face loop
261
262
263 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
264 // Step 3: ELEMENTS AND ELEMENT MATRICES ON BOUNDARY CORNERS
265 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
266 for (unsigned int corner = 0; corner < n_corner.getNx() ;corner++ ) {
267 element = n_corner[corner];
268 cellType = "corner";
269 cell.classifyNodes(dofNode,&nDofNodes,&nBoundaryNodes,cellType,corner);
270 unsigned int ix=0;
271 unsigned int iy=0;
272 unsigned int iz=0;
273
274 if (corner<=3){
275 iz=0;
276 if (corner==0 || corner==3){
277ix=0;
278 }else{
279ix=nx-2;
280 }
281 if (corner==2 || corner==3){
282iy=ny-2;
283 } else{
284iy=0;
285 }
286 }else if (corner<=7){
287 iz=nz-1;
288 if (corner==4 || corner==7){
289ix=0;
290 }else{
291ix=nx-2;
292 }
293 if (corner==6 || corner==7){
294iy=ny-2;
295 } else{
296iy=0;
297 }
298 }
299
300 computeDriftVector(b, ix, iy, iz);
301
302 computeDiffusionTensor(K, ix, iy, iz);
303
304 computeElementParameters(element,DETERMINANTJ,J0M,J0,J,b,K,deltak,qualla,ix,iy,iz);
305
306 //DirichletBoundaryValues hasn't been implemented
307 computeDirichletBoundaryValues( element,DETERMINANTJ,J0M,J0,b,K,deltak,dofNode,nDofNodes,nBoundaryNodes,var00,Psi_loc);
308 }//end corner loop
309
310 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
311 // Step 4: ELEMENTS AND ELEMENT MATRICES ON BARS
312 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
313 size_t nz_bar_size, tmp2;
314 nz_bar.size(nz_bar_size, tmp2);
315 //for (int bar = 0; bar < nz_bar.size(); bar++ )
316 for (size_t bar = 0; bar < nz_bar_size; bar++ ){
317for (int layer=nz_bar[bar][0];layer<=nz_bar[bar][1];layer++){
318 for (int depth=ny_bar[bar][0];depth<=ny_bar[bar][1];depth++) {
319for (int width=nx_bar[bar][0];width<=nx_bar[bar][1];width++){
320
321element = layer*(nx-1)*(ny-1) + depth*(nx-1) + width;
322computeDriftVector(b, width, depth, layer);
323computeDiffusionTensor(K,width, depth, layer);
324cellType = "bar";
325cell.classifyNodes(dofNode,&nDofNodes,&nBoundaryNodes,cellType,(int)bar);
326
327computeElementParameters(element, DETERMINANTJ,J0M,J0,J,b,K,deltak,qualla,width,depth,layer);
328//DirichletBoundaryValues hasn't been implemented >>remove?
329computeDirichletBoundaryValues( element,DETERMINANTJ,J0M,J0,b,K,deltak,dofNode,nDofNodes,nBoundaryNodes,var00,Psi_loc);
330}// end width loop
331 }// end depth loop
332}// end layer loop
333 } // end bar loop
334}
335
336/**
337 * @brief Compute diffusion tensor
338 * Computes the diffusion tensor of an element, accesses the nodes array and the nodeMap array
339 * @param element element number
340 * @param K diffusion coefficient
341 * @param layer layer of the element
342*/
343void SnowDriftA3D::computeDiffusionTensor(double K[3][3], const unsigned int ix, const unsigned int iy, const unsigned int iz)
344{
345 double Length_Scale;
346 double K_model,z_loc,K_min,du;
347 double K_max=1.6;
348
349 int version=1; //1: old version with mixing length, 2: take K from ARPS, 3: version where K depends on standard deviation of wind speed within a layer, 4: K is calculated after the windfield is read (see SnowDrift.cc),calulation was moved because it is time consuming and is now calculated only once.
350 if (READK){
351 version=2;
352 }
353
354 for (int i = 0; i < 3; i++ ) {
355
356for (int j = 0; j < 3; j++ ){
357 K[i][j] = 0;
358}
359 }
360
361 K_model = 0.; K_min = 0.;
362
363 if (version==1){
364//loop over 4 nodes
365unsigned int ii, jj;
366for ( ii=ix; ii<=ix+1; ii++){
367 for (jj = iy; jj<=iy+1; jj++) {
368
369Length_Scale = fabs(nodes_z.grid3D(ii,jj,iz)-nodes_z.grid3D(ii,jj,iz+1));
370
371//vertical height
372z_loc = (nodes_z.grid3D(ii,jj,iz) - nodes_z.grid3D(ii,jj,0)) + 0.5*Length_Scale;
373
374/* Approximation of ABL neutral: 1/8*dudz*SQR(1/((1/k(z+z0))+(1/LengthScale))) */
375du = sqrt( Optim::pow2( nodes_u.grid3D(ii,jj,iz+1) - nodes_u.grid3D(ii,jj,iz))
376 + Optim::pow2( nodes_v.grid3D(ii,jj,iz+1) - nodes_v.grid3D(ii,jj,iz) )
377 + Optim::pow2( nodes_w.grid3D(ii,jj,iz+1) - nodes_w.grid3D(ii,jj,iz) )
378 );
379
380K_model += 0.25*du/Length_Scale*(Optim::pow2(1./((1./(0.4*z_loc))+(1./Length_Scale))));
381
382/* Approximation of k*ustar*z for the surface layer */
383K_min += 0.25*0.4*z_loc*0.1;
384 }
385}
386 K_min = std::max(0.1,K_min);
387 }else if (version==2){
388 //use the diffusivity as given in ARPS >> will only work if you use wind fields of ARPS with a somewhat longer integration step (thus incl. turbulence). Some tests with K of mean //flow fields > solver couldn't find steady state of c
389 double temp_K=0.;
390 unsigned int ii,jj;
391 for ( ii=ix; ii<=ix+1; ii++){
392for (jj = iy; jj<=iy+1; jj++) {
393 temp_K=temp_K+nodes_K.grid3D(ii,jj,iz)+nodes_K.grid3D(ii,jj,iz+1);
394
395}
396 }
397 //now we have added all possible K of the 8 nodes around element, take the average:
398 K_model=temp_K/8.;
399 K_min=0.1;
400
401
402 }
403 //make sure K_model is in range min-max
404 K_model=std::max(K_min,K_model);
405 K_model=std::min(K_max,K_model);
406
407 K[2][2] = K_model;
408 K[1][1] = K_model;
409 K[0][0] = K_model;
410 //extra test
411 K[0][0]=std::min(std::max(K[0][0],K_min),50.*K_model);
412 K[1][1]=std::min(std::max(K[1][1],K_min),50.*K_model);
413 K[2][2]=std::min(std::max(K[2][2],K_min),10.*K_model);//ori 10
414}
415
416/**
417 * @brief Compute the drift vector of an element
418 * @param element element number
419 * @param b drift vector
420 */
421void SnowDriftA3D::computeDriftVector(double b[3], const unsigned int ix, const unsigned int iy, const unsigned int iz)
422{
423 //reset
424 for (int i = 0; i < 3; i++) {
425 b[i]=0.;
426 }
427
428 unsigned int ii, jj, kk;
429 for ( ii=ix; ii<=ix+1; ii++){
430for (jj = iy; jj<=iy+1; jj++) {
431 for (kk = iz; kk<=iz+1; kk++) {
432b[0]+=1/8.*nodes_u.grid3D(ii,jj,kk);
433b[1]+=1/8.*nodes_v.grid3D(ii,jj,kk);
434b[2]+=1/8.*(nodes_w.grid3D(ii,jj,kk)-nodes_wstar.grid3D(ii,jj,kk));
435 }
436}
437 }
438}
439
440/**
441 * @brief Compute parameters of the SUPG method
442 * Compute the parameters of the SUPG method for each element
443 * @param element Element number
444 * @param DETERMINANTJ Determinant of Jacobian matrix
445 * @param J0M
446 * @param J0 JO matrix
447 * @param J Jacobian matrix
448 * @param b drift vector
449 * @param K Diffusion tensor
450 * @param &deltak parameter SUPG method
451 * @param &qualla parameter SUPG method (to vary delta_k)
452 */
453void SnowDriftA3D::computeElementParameters(const int& element,
454 double DETERMINANTJ[8],
455 double J0M[3][3][8],
456 double J0[3][3],// the J0 matrix
457 double J[3][3],// the Jacobian matrix
458 double b[3],
459 double K[3][3],
460 double &deltak,
461 double &qualla_loc,
462 const int ix,
463 const int iy,
464 const int iz)
465{
466 //cad is a dummy vector the integration point for the Gaussian
467 //quadrature is assigned to
468 double cad[3];
469 double epsilon;
470 double hk; // the diameter (=longest side) of the respective
471 // element
472 double Pe; // the local Peclet number
473
474 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
475 //element infos
476 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
477 epsilon = 1/3.* (K[0][0]+K[1][1]+K[2][2]);
478
479 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
480 //hk
481 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
482 hk=fabs(nodes_x.grid3D(ix,iy,iz+1)-nodes_x.grid3D(ix+1,iy,iz+1));
483
484 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
485 // deltak
486 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
487 Pe = sqrt(b[0]*b[0]+b[1]*b[1]+b[2]*b[2]) * hk/12/epsilon; /*calculate the local Peclet number*/
488 if (Pe<1){
489 deltak=hk*hk/12/epsilon*qualla_loc; //qualla is parameter to vary delta_k, standard:1
490 }else{
491 deltak=hk/2/sqrt(b[0]*b[0]+b[1]*b[1]+b[2]*b[2])*qualla_loc;
492 }
493
494 for (int i = 0 ; i < 8; i++){ //loop over all 8 points
495 cad[0] = qPoint(0, i);
496 cad[1] = qPoint(1, i);
497 cad[2] = qPoint(2, i);
498
499 // this function will calculate the Jacobian of the
500 // isoparametric point transformation on the
501 // respective element and stores its determinant in
502 // DETERMINANTJ[i] for each of the integration
503 // points and store the Jacobian J at the respective
504 // point in J; //note that you'll only use J once
505 // more for J0
506 Jacobian(DETERMINANTJ,J,element,cad,i,ix,iy,iz);
507
508
509 // this function will calculate the J0 matrix at the
510 // integration point and store it in the variable
511 // J0*/
512 J0fun(J0,J);
513
514 for (int k = 0; k < 3; k++ ){
515 for (int j = 0; j < 3; j++){
516 J0M[k][j][i] = J0[k][j];
517 }
518 }
519 }
520}
521
522
523/**
524 * @brief compute element system
525 * ....
526 * @param element element number
527 * @param nDofNodes degrees of freedom of element
528 * @param Ael
529 * @param Del
530 * @param stationary
531 * @param DETERMINANTJ
532 * @param J0M
533 * @param b wind field
534 * @param K diffusion tensor
535 * @param deltak parameter
536 * @param dt
537 * @param f source
538 * @param Psi
539*/
540void SnowDriftA3D::computeElementSystem(int &element,
541 int &nDofNodes,
542 int* dofNode,
543 double Ael[9][9],
544 double Del[9][9],
545 bool stationary,
546 double DETERMINANTJ[8],
547 double J0M[3][3][8],
548 double b[3],
549 double K[3][3],
550 double &deltak,
551 const double &dt,
552 CDoubleArray& f_loc,
553 CDoubleArray& Psi_loc)
554{
555 //element matrices
556 double Bel[9][9];
557 double Cel[9][9];
558 double Apdxel[9][9];
559 double Adxdxel[9][9];
560
561 //reset element matrices
562 for (int i = 0; i<9; i++ ){
563 for (int j = 0; j<9; j++ ){
564Ael[i][j] = 0; //LH
565Bel[i][j] = 0;
566Cel[i][j] = 0;
567Apdxel[i][j] = 0;
568Adxdxel[i][j] = 0;
569 }
570 }
571
572 //------------------------------------
573 // calculation of the element matrices
574 //------------------------------------
575 for (int i = 0; i < nDofNodes ; i++ ){
576 double auxSum = 0;
577 for (int j = 0; j < nDofNodes ; j++){
578//element matrices are indexed from 1..8
579Bel[ 1+dofNode[i] ][ 1+dofNode[i] ] = GQIntB(DETERMINANTJ,dofNode[i],dofNode[j]);
580Cel[ 1+dofNode[i] ][ 1+dofNode[j] ] = GQIntC(DETERMINANTJ,J0M,dofNode[i],dofNode[j],b,K);
581Adxdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] = GQIntAdxdx(DETERMINANTJ,J0M,dofNode[i],dofNode[j],b,deltak);
582Apdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] = GQIntApdx(DETERMINANTJ,J0M,dofNode[i],dofNode[j],b,deltak);
583
584//assemble matrix according to crank nicolson
585if ( stationary ) {//LH
586Ael[ 1+dofNode[i] ][ 1+dofNode[j] ] =
587( Adxdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] + Cel[ 1+dofNode[i] ][ 1+dofNode[j] ] );
588}else{
589Ael[ 1+dofNode[i] ][ 1+dofNode[j] ] =
590( Apdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] + Bel[ 1+dofNode[i] ][ 1+dofNode[j] ] ) / dt
591+ ( Adxdxel[ 1+dofNode[i] ][ 1+dofNode[j]] + Cel[ 1+dofNode[i] ][ 1+dofNode[j] ] ) * 0.5;
592}
593
594//assemble right hand matrix according to crank nicolson
595if ( stationary ){//LH
596Del[ 1+dofNode[i] ][ 1+dofNode[j] ] = 0;
597} else {
598 Del[ 1+dofNode[i] ][ 1+dofNode[j] ] =
599( Apdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] + Bel[ 1+dofNode[i] ][ 1+dofNode[j] ] ) / dt
600- ( Adxdxel[ 1+dofNode[i] ][ 1+dofNode[j]] + Cel[ 1+dofNode[i] ][ 1+dofNode[j] ] ) * 0.5;
601}
602
603//assemble right hand side
604// if ( stationary )//LH
605// {
606 auxSum += ( Bel[ 1+dofNode[i] ][ 1+dofNode[j] ] + Apdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] )
607 * f_loc[ nodeMap[element][ dofNode[j] ] ] ;
608
609// }
610// else
611// {
612// auxSum += ( ( ( Apdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] + Bel[ 1+dofNode[i] ][ 1+dofNode[j] ] )/dt
613// - ( Cel[ 1+dofNode[i] ][ 1+dofNode[j] ] + Adxdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] )*0.5
614// ) * c[ nodeMap[element][ dofNode[j] ] ]
615// + ( Bel[ 1+dofNode[i] ][ 1+dofNode[j] ] + Apdxel[ 1+dofNode[i] ][ 1+dofNode[j] ] )
616// * f[ nodeMap[element][ dofNode[j] ] ]
617// );
618// }
619 }
620 Psi_loc[ nodeMap[element][ dofNode[i] ] ] += auxSum;
621 }
622}
623
624/**
625 * @brief Add element matrix
626 * given a sparse matrix A in CSR format ( specified by
627 * sA, colInd, rowPtr) and the element matrix Bel, this function adds
628 * the element contributions to the global ones ; the first
629 * length_spec entries of the vector spec contains the indices of the
630 * matrix which are inserted
631 * Comments : Beware of the enumeration of the matrix sA (size is one
632 * bigger than required, zero index is unused => dangerous
633 * accesses nodeMap, nodeMap
634 * @param sA
635 * @param colInd
636 * @param rowPtr
637 * @param Bel
638 * @param element
639 * @param spec
640 * @param length_spec
641*/
642void SnowDriftA3D::addElementMatrix( CDoubleArray& sA_loc,
643 const CIntArray& colInd,
644 const CIntArray& rowPtr,
645 const double Bel[9][9],
646 const int element,
647 const int* spec,
648 const int length_spec)
649{
650 //global indices
651 int I;
652 int J;
653 for ( int i = 0; i < length_spec; i++){
654 I = nodeMap[element][ spec[i] ];
655
656 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
657 //insert entry (i,i)
658 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
659 //binary search of I's column index
660 int first = rowPtr[I];
661 int last = rowPtr[I+1] - 1;
662 int mid;
663 while ( first <= last ){
664 mid = (first + last) / 2; // compute mid point.
665 if ( I > colInd[mid] ){
666first = mid + 1; // repeat search in top half.
667 }else if ( I < colInd[mid] ) {
668last = mid - 1; // repeat search in bottom half.
669 }else {
670sA_loc[ mid ] += Bel[ spec[i]+1 ][ spec[i]+1 ]; // got it
671 break;
672 }
673 }
674
675 for ( int j = i + 1; j < length_spec; j++){
676 J = nodeMap[element][ spec[j] ];
677 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
678 //insert entry (i,j)
679 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
680
681 //binary search of J's column index
682 first = rowPtr[I];
683 last = rowPtr[I+1] - 1;
684 while (first <= last){
685mid = (first + last) / 2; // compute mid point.
686if ( J > colInd[mid] ){
687first = mid + 1; // repeat search in top half.
688}else if ( J < colInd[mid] ) {
689last = mid - 1; // repeat search in bottom half.
690}else{
691sA_loc[ mid ] += Bel[ spec[i]+1 ][ spec[j]+1 ]; // got it
692break;
693}
694 }
695
696 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
697 //insert entry (j,i)
698 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
699 //binary search of I's column index
700 first = rowPtr[J];
701 last = rowPtr[J+1] - 1;
702 while ( first <= last ) {
703 mid = (first + last) / 2; // compute mid point.
704 if ( I > colInd[mid] ) {
705 first = mid + 1; // repeat search in top half.
706 }else if ( I < colInd[mid] ) {
707 last = mid - 1; // repeat search in bottom half.
708 }else {
709 sA_loc[ mid ] += Bel[ spec[j]+1 ][ spec[i]+1 ]; // got it
710 break;
711 }
712 }
713}//i
714 }//j
715}
716
717//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
718// Function : computeDirichletBoundaryValues
719// Authors : Marc Ryser
720//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
721void SnowDriftA3D::computeDirichletBoundaryValues(int element,
722 double DETERMINANTJ[8],
723 double J0M[3][3][8],
724 double J0[3][3],// the J0 matrix
725 double b[3],
726 double K[3][3],
727 double deltak,
728 int node[8],
729 int nDofNodes,
730 int nBoundaryNodes,
731 CDoubleArray& c00_loc,
732 CDoubleArray& Psi_loc )
733{
734//since this method is currently empty, we just remove the warnings...
735//TODO: implement this method or get rid of it
736(void)element;
737(void)DETERMINANTJ;
738(void)J0M;
739(void)J0;
740(void)b;
741(void)K;
742(void)deltak;
743(void)node;
744(void)nDofNodes;
745(void)nBoundaryNodes;
746(void)c00_loc;
747(void)Psi_loc;
748 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
749 // incorporation of inhomogeneous boundary conditions
750 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
751
752// for (int i = 0; i < nDofNodes; i++ )
753// {
754// for (int j = nDofNodes; j < nDofNodes + nBoundaryNodes; j++)
755// {
756// Psi[ nodeMap[ element][ node[i] ] ] +=
757// c00[ nodeMap[ element][ node[j] ] ]
758// * GQIntC( DETERMINANTJ,J0M, node[i], node[j],b,K);
759
760// Psi[ nodeMap[element][ node[i] ] ] +=
761// c00[ nodeMap[element][ node[j] ] ]
762// * GQIntAdxdx(DETERMINANTJ,J0M, node[i], node[j],b,deltak);
763
764// rhs[ nodeMap[element][ node[j] ] ] = c00[ nodeMap[element][ node[j] ] ];
765// // printf("c00[%d]=%lf\n", nodeMap[element][ node[j]], c00[ nodeMap[element][ node[j] ] ] );
766// zeroRow( nodeMap[element][ node[j] ] );
767// }
768// }
769
770 //if bottom
771// if ( element <= (nx-1)*(ny-1) )
772// {
773 //enumeration from 1..8
774// double BCel[9][9];
775// double surfaceMetric;
776// double qp[3];
777// double PHI[8];
778
779// for ( int i = 0; i < 8; i++)
780// {
781// for ( int j = 0; j < 8; j++)
782// {
783// BCel[i+1][j+1] = 0;
784
785// //determine element normal to bottom
786// for ( int k = 0; k < 4; k++ )
787// {
788// }
789
790// //loop over quadrature points
791// for ( int k = 0; k < 4; k++ )
792// {
793// qp[0]= qPoint[0][k];
794// qp[1]= qPoint[1][k];
795// qp[2]= -1;
796
797// phi(PHI,qp);
798// Jacobian(DETERMINANTJ,J,element,qp,k);
799// J0fun(J0,J);
800// surfaceMetric = sqrt( J0[1][3]*J0[1][3]
801// +J0[2][3]*J0[2][3]
802// +J0[3][3]*J0[3][3] );
803
804
805
806
807
808// BCel[i+1][j+1] += ( surfaceMetric
809// * gamma
810// * b[2]
811// * (PHI[i]-gDirichlet[ nodeMap[element][k]] )
812// + gNeumann[ nodeMap[element][k] ]
813// ) * PHI[j];
814
815// }
816// }
817// }
818// addElementMatrix(sA,colA,rowA,BCel,element,dofNode,nDofNodes+nBoundaryNodes);//LH
819// // }
820
821
822
823}
824
825
826
827/**
828 * @brief compute deposition flux
829 * Authors : Marc Ryser
830 * Description : calculates the deposition flux i.e. (b c + K grad c) at the
831 * artificial layer of nodes and write it to the flux_x,y,z variable
832 * Comments : Beware of the enumeration of elements and nodes, Must be
833 * consistent with SnowDriftA3D::SnowMassChange(...) Note: this solution
834 * has temporary character since flux_x has size (nx-1)*(ny-1) and
835 * therefore only internal fluxes are written
836 * @param CDoubleArray&c Concentration of snow
837 * @param theta height above bottom of element
838*/
839void SnowDriftA3D::computeDepositionFlux(const CDoubleArray& concentration,
840 const double theta_l)
841{
842 double DETERMINANTJ[8];
843 double J0[3][3];
844 double J[3][3];
845 double TT[3][8];
846// double b[3];
847 double K[3][3];
848
849 //cad is a dummy vector the integration point for the Gaussian
850 //quadrature is assigned to
851 double cad[3];
852
853 double hivec[3];
854 double gradc[3];
855
856 int element;
857
858 // loop over all the second layer
859 // elements having the bottom node 0
860 // in the interior of the domain
861 for (unsigned int k = 1; k < ny-1; k++ ){
862 for (unsigned int l = 1; l < nx-1; l++){
863 // the respective element of the second layer
864 element = (nx-1)*(ny-1) + k*(nx-1) + l;
865// int element0 = k*(nx-1) + l;
866 int node = nodeMap[element][0]-nx*ny;
867 //reset fluxes
868 flux_x[node] = 0;
869 flux_y[node] = 0;
870 flux_z[node] = 0;
871
872 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
873 // advective contribution
874 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
875flux_x[node] += concentration[nodeMap[element][0]] * nodes_u.grid3D(l,k,1);
876flux_y[node] += concentration[nodeMap[element][0]] * nodes_v.grid3D(l,k,1);
877flux_z[node] += concentration[nodeMap[element][0]] * ( nodes_w.grid3D(l,k,1)-nodes_wstar.grid3D(l,k,1));
878
879 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
880 // diffusive contribution
881 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
882 computeDiffusionTensor(K, l,k,1); //take layer consistent with element
883
884 // This Is The Vector Containing The Point Where We Want To
885 // Calculate K gradc dot n theta is the height above bottom of
886 // the element, very small if possible to get as close as
887 // possible to the point
888 cad[0] = 1 - theta_l;
889 cad[1] = -1 + theta_l;
890 cad[2] = -1 + theta_l;
891
892 // calculation of J at the point cad and storing the determinant in the
893 // 0th (=kth component) component of DETERMINANTJ
894 Jacobian( DETERMINANTJ, J, element, cad, 0 ,l,k,1);
895
896 // calculation of J0 at the point cad
897 J0fun(J0,J);
898
899 //calculation of the gradient of all 8 shape functions at the point cad in the reference element
900 TTfun(TT,cad);
901
902 for (int j = 0; j < 3; j++){ // the 3 components of gradc
903 gradc[j] = 0;
904 for (int m = 0; m < 8; m++){ //add over all 8 shape functions
905hivec[j] = 0;// need to clear hivec[j] for the matrix product of J0 with gradphi(m)
906for (int n = 0; n < 3; n++){ // dummy loop for matrix product
907 hivec[j] += J0[j][n] * TT[n][m];//store J0 * gradphi(m)
908}
909gradc[j] += concentration[nodeMap[element][m]] * (1/DETERMINANTJ[0]) * hivec[j];//gradc=sum_m
910//(c(m)*1/detJ
911//*
912//J0*gradphi(m))
913 }
914 }
915 // k=1, l=2, and element given above implies
916 // nodeMap[element][0]=1, ie flux_x is filled up from the index 1
917 // k=ny-2, l=nx-1, and element given above implies
918 // nodeMap[element][0]=3721=(nx-2)*(ny-2)
919 // note: flux_x,y,z have size (nx-1)*(ny-1) such that the
920 // remaining memory is unaltered (better: resize flux_x)
921 flux_x[ node ] -= K[0][0] * gradc[0];
922 flux_y[ node ] -= K[1][1] * gradc[1];
923 flux_z[ node ] -= K[2][2] * gradc[2];
924} //end of kicking all elements on 1st interior element layer in l!!!
925 } //dito in k!!
926}
927
928/**
929 * @brief compute deposition flux
930 * Authors : Marc Ryser
931 * Description : calculates the deposition flux i.e. (b c + K grad c) at the
932 * artificial layer of nodes and write it to the flux_x,y,z variable
933 * Comments : Beware of the enumeration of elements and nodes, Must be
934 * consistent with SnowDriftA3D::SnowMassChange(...) Note: this solution
935 * has temporary character since flux_x has size (nx-1)*(ny-1) and
936 * therefore only internal fluxes are written
937 * @param CDoubleArray&c Concentration of snow
938 * @param theta height above bottom of element
939*/
940void SnowDriftA3D::computeDepositionFluxSublimation(const CDoubleArray& concentration,
941 const double theta_l)
942{
943 double DETERMINANTJ[8];
944 double J0[3][3];
945 double J[3][3];
946 double TT[3][8];
947// double b[3];
948 double K[3][3];
949
950 //cad is a dummy vector the integration point for the Gaussian
951 //quadrature is assigned to
952 double cad[3];
953
954 double hivec[3];
955 double gradc[3];
956
957 int element;
958
959 // loop over all the second layer
960 // elements having the bottom node 0
961 // in the interior of the domain
962 for (unsigned int k = 1; k < ny-1; k++ ){
963 for (unsigned int l = 1; l < nx-1; l++){
964 // the respective element of the second layer
965 element = (nx-1)*(ny-1) + k*(nx-1) + l;
966// int element0 = k*(nx-1) + l;
967 int node = nodeMap[element][0]-nx*ny;
968 //reset fluxes
969 flux_x_subl[node] = 0;
970 flux_y_subl[node] = 0;
971 flux_z_subl[node] = 0;
972
973 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
974 // advective contribution
975 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
976 flux_x_subl[node] += concentration[nodeMap[element][0]] * nodes_u.grid3D(l,k,1);
977 flux_y_subl[node] += concentration[nodeMap[element][0]] * nodes_v.grid3D(l,k,1);
978 flux_z_subl[node] += concentration[nodeMap[element][0]] * ( nodes_w.grid3D(l,k,1)-nodes_wstar.grid3D(l,k,1));
979
980 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
981 // diffusive contribution
982 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
983 computeDiffusionTensor(K, l,k,1); //take layer consistent with element
984
985 // This Is The Vector Containing The Point Where We Want To
986 // Calculate K gradc dot n theta is the height above bottom of
987 // the element, very small if possible to get as close as
988 // possible to the point
989 cad[0] = 1 - theta_l;
990 cad[1] = -1 + theta_l;
991 cad[2] = -1 + theta_l;
992
993 // calculation of J at the point cad and storing the determinant in the
994 // 0th (=kth component) component of DETERMINANTJ
995 Jacobian( DETERMINANTJ, J, element, cad, 0 ,l,k,1);
996
997 // calculation of J0 at the point cad
998 J0fun(J0,J);
999
1000 //calculation of the gradient of all 8 shape functions at the point cad in the reference element
1001 TTfun(TT,cad);
1002
1003 for (int j = 0; j < 3; j++){ // the 3 components of gradc
1004 gradc[j] = 0;
1005 for (int m = 0; m < 8; m++){ //add over all 8 shape functions
1006hivec[j] = 0;// need to clear hivec[j] for the matrix product of J0 with gradphi(m)
1007for (int n = 0; n < 3; n++){ // dummy loop for matrix product
1008 hivec[j] += J0[j][n] * TT[n][m];//store J0 * gradphi(m)
1009}
1010gradc[j] += concentration[nodeMap[element][m]] * (1/DETERMINANTJ[0]) * hivec[j];//gradc=sum_m
1011//(c(m)*1/detJ
1012//*
1013//J0*gradphi(m))
1014 }
1015 }
1016 // k=1, l=2, and element given above implies
1017 // nodeMap[element][0]=1, ie flux_x is filled up from the index 1
1018 // k=ny-2, l=nx-1, and element given above implies
1019 // nodeMap[element][0]=3721=(nx-2)*(ny-2)
1020 // note: flux_x,y,z have size (nx-1)*(ny-1) such that the
1021 // remaining memory is unaltered (better: resize flux_x)
1022 flux_x_subl[ node ] -= K[0][0] * gradc[0];
1023 flux_y_subl[ node ] -= K[1][1] * gradc[1];
1024 flux_z_subl[ node ] -= K[2][2] * gradc[2];
1025} //end of kicking all elements on 1st interior element layer in l!!!
1026 } //dito in k!!
1027}
1028
1029void SnowDriftA3D::zeroRow(int node)
1030{
1031 //version 3
1032 int I = node;
1033 int nnz = rowA[I+1]-rowA[I];
1034
1035 //loop over all indices in row I
1036 for (int k = 0; k < nnz; k++){
1037 int J = colA[ rowA[I] + k ];
1038 //set diagonal element to 1
1039 if ( I == J ){
1040 sA[rowA[I] + k] = 1;
1041 }else{
1042 //zero off diagonal elements
1043 //IJ
1044 sA[rowA[I] + k] = 0;
1045
1046 //JI
1047 //binary search of I's column index
1048 int first = rowA[J];
1049 int last = rowA[J+1] - 1;
1050 int mid;
1051 while (first <= last) {
1052 mid = (first + last) / 2; // compute mid point.
1053 if ( I > colA[mid] ) {
1054 first = mid + 1; // repeat search in top half.
1055 } else if ( I < colA[mid] ) {
1056 last = mid - 1; // repeat search in bottom half.
1057 }else{
1058 sA[ mid ] = 0; // got it
1059 break;
1060 }
1061 }
1062 }
1063 }
1064}

Archive Download this file

Revision: HEAD