Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/snowdrift/Suspension.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/*------------------------------------------------------------------------------------------+
19 | This module is the FINITE ELEMENT solution to Snow Diffusion |
20 +------------------------------------------------------------------------------------------*/
21/********************************************************************************************/
22#include <alpine3d/snowdrift/SnowDrift.h>
23#include <alpine3d/AlpineMain.h>
24#include <iostream>
25
26using namespace mio;
27using namespace std;
28
29#define DRIFT_MARC 1
30
31void SnowDriftA3D::Diffusion(double deltaT,double &diff_max, double t1)
32{//OLD VERSION
33(void)deltaT;
34(void)diff_max;
35(void)t1;
36std::cerr<<"[E] This method is obsolete, it should not be used!!"<<std::endl;
37}
38 /* End of Diffusion */
39
40void SnowDriftA3D::Suspension()
41{
42 const double eps=0.00005;
43 double diff_max[2];
44 diff_max[1] = 0.0;
45 double t_diff=0.0;
46 static int tmpindex=0;
47
48 if (DRIFT_MARC)
49 {
50int timeStep;
51const double maxTime = 0.8; // the time and upper limit of the time interval of integration
52const int maxTimeStep = (int) ceil( maxTime/dt_diff );
53
54// Inititalize system
55initializeSystem( colA, rowA, sA,sB,rhs,f,Psi,c,c00, CON);
56prepareSolve();
57timeStep=0;
58SolveEquation(timeStep, maxTimeStep,CON);
59values_elements_to_nodes(nodes_c,c);
60double nodes_min=nodes_c.grid3D.getMin();
61if (nodes_min<0){
62cout<<"!!!WARNING!!! at least one point with negative snow concentration! (minimum c="<<nodes_min<<")"<<endl;
63}
64
65// compute deposition flux
66computeDepositionFlux( c, theta );
67
68if (SUBLIMATION && (nodes_c.grid3D.getMean()>0.00000)){
69std::cout<<"Start sublimation routine"<< std::endl;
70
71//now find the steady state sublimation
72iterativeSublimationCalculation(timeStep, maxTimeStep);
73
74//steady state sublimation is found, if the concentration feedback is switched off, the new concentration field will have to be calculated here to be able to consider the effect of sublimation on concentration and deposition
75if (C_FB== 0){
76 //recalculate concentration field
77 initializeSystem( colA, rowA,sA,sB, rhs, f, Psi, c,c00, CON);
78 //change source f to steady state sublimation
79 Grid3DObject tmpSink_c;
80 tmpSink_c=nodes_Subl_ini;
81 tmpSink_c.grid3D=nodes_Subl_ini.grid3D*(-1);
82 values_nodes_to_elements(tmpSink_c,f);
83
84 prepareSolve();
85 SolveEquation(0, maxTimeStep, CON);
86 values_elements_to_nodes(nodes_c,c);
87 nodes_min=nodes_c.grid3D.getMin();
88 if (nodes_min<0){
89cout<<"!!WARNING!! at least one point with negative concentration! (minimum c="<<nodes_min<<")"<<endl;
90 }
91}
92computeDepositionFluxSublimation( c, theta );
93}
94 }//driftmarc
95 else {
96do {
97t_diff += dt_diff;
98diff_max[0] = diff_max[1];
99/* Calculate the diffusion */
100Diffusion(dt_diff, diff_max[1], t_diff);
101printf("\n Local c: %f Diff Time: %f s diff_max:%f\n", nodes_c.grid3D(2,2,2),t_diff,diff_max[1]);
102
103tmpindex++;
104//DEBUG("SEQ FLUX (ix=31): x=%f, y=%f, z=%f",
105//checksum(flux_x,31, nx-1),checksum(flux_y,31, nx-1),checksum(flux_z,31, nx-1));
106//DEBUG("SEQ FLUX (ix=33): x=%f, y=%f, z=%f",
107//checksum(flux_x,33, nx-1),checksum(flux_y,33, nx-1),checksum(flux_z,33, nx-1));
108
109} while ( (diff_max[1] > eps) && (t_diff/dt_diff < 100.0) );
110
111printf("\n Converged after time %f\n", t_diff);
112}
113//DEBUG("DRIFT FLUX #%d: x=%f, y=%f, z=%f", tmpindex, checksum(flux_x), checksum(flux_y), checksum(flux_z));
114}
115
116
117/**
118 * @brief Solve the advection diffusion equation
119 * @param timeStep (parameter only for non stationary)
120 * @param maxTimeStep (parameter only for non stationary)
121 * @param param variable that should be solved
122 */
123void SnowDriftA3D::SolveEquation(int timeStep, int maxTimeStep,const param_type param)
124{
125// Solve equation
126if ( !(STATIONARY) ) {
127
128 while ( timeStep < maxTimeStep ) {
129if ( timeStep == 0 ) {
130 //extract the initial conditions for time=0!
131 for (unsigned int kk = 0; kk<= nz-3; kk++) {
132// kick all the internal node layers!
133for (unsigned int jj = 0; jj <= ny-3; jj++ ) {
134 for (unsigned int ii = 1; ii <= nx-2; ii++ ) {
135c[kk*(nx-2)*(ny-2)+jj*(nx-2)+ii] = c00[(kk+1)*nx*ny+(jj+1)*nx+ii];
136 }
137}
138 }
139}
140
141if (param==HUM){
142matmult(rhs,q,sB,colA,rowA);
143} else if (param==CON){
144matmult(rhs,c,sB,colA,rowA);
145} else if (param==TEM){
146matmult(rhs,T,sB,colA,rowA);
147}
148
149//construct the rhs of the system to solve: Ac=rhs
150for (unsigned int i = 0; i < nDOF; i++) {
151rhs[i] += ( -Psi[i] );
152}
153
154//solve system
155double test;
156if (param==HUM){
157bicgStab(q,rhs,sA,colA,rowA,900,1e-5, test);
158} else if (param==CON){
159bicgStab(c,rhs,sA,colA,rowA,900,1e-5, test);
160} else if (param==TEM){
161bicgStab(T,rhs,sA,colA,rowA,900,1e-5, test);
162}
163//update time
164timeStep++;
165 }
166} else {
167//construct the right hand side rhs of the system A*c = rhs
168for (unsigned int i = 0; i < nDOF; i++ ) {
169rhs[i] += ( -Psi[i] );
170
171//LH_BC
172// rhs[i] = 0.0005;
173//rhs[i] = res1[i] - Psi[i] + res2[i];//Here real line
174}
175
176double testres=10.;
177CDoubleArray saveRhs=rhs;
178CIntArray saveRowA=rowA;
179if (param==HUM){
180 std::cout<<"Solving equation for q with rhs=" <<rhs[18]<<std::endl;
181 bicgStab(q,rhs,sA,colA,rowA,3000,1e-11,testres);//call the biconjgrad method,ori 1e-10
182} else if (param==CON){
183 std::cout<<"Solving equation for c with rhs=" <<rhs[18]<<std::endl;
184 /*std::string fname = std::string("../results/") + calcDate.toString(Date::NUM)+ "0" + std::string(".c00");
185writefield(fname,HUM);*/
186 bicgStab(c,rhs,sA,colA,rowA,3000,1e-10,testres); //ori 1e-10, 1000
187} else if (param==TEM){
188 std::cout<<"Solving equation for T with rhs=" <<rhs[18]<<std::endl;
189 bicgStab(T,rhs,sA,colA,rowA,3000,1e-11,testres);
190}
191}
192}
193
194/**
195@brief
196Calculate the steady state sublimation in several steps.
197Single feedbacks can be switched on or off in SnowDrift.h
198*/
199void SnowDriftA3D::iterativeSublimationCalculation(int timeStep, int maxTimeStep)
200{
201double factorf=(-1); //fraction of initial sublimation that should be used to solve the equation. -1 by definition
202
203//Initial field of concentration has been calculated, now calculate the initial sublimation
204Sublimation();
205
206//initialize some parameters
207int count=0; //count iterations
208double maxdif_s=0.; //the maximum difference in sublimation between 2 iterations
209double maxdif_c=0.; //the maximum difference in concentration between 2 iterations
210double test=0; //controls whether iteration stops
211double mean_startsubl=0.;//mean sublimation at start iteration
212double mean_endsubl=0; //mean final sublimation
213
214//initialize init sublimation
215nodes_Subl_ini.grid3D=0.;
216
217//Save initial sublimation and adjust to start iteration with a value probably closer to steady-state sublimation
218if (C_FB==1 || Q_FB==1 || T_FB==1){
219 nodes_Subl_ini.grid3D=nodes_Subl.grid3D*0.8;
220}else{
221 nodes_Subl_ini=nodes_Subl;
222}
223mean_startsubl=nodes_Subl.grid3D.getMean();
224
225//now start the iteration
226while (count<8 && test==0 && mean_startsubl!=0.){ //Usually not more than 7 iterations are needed but can be adjusted. Test is defined at the end of the loop and will test changes in sublimation and concentration
227count=count+1;
228std::cout <<"Sublimation iteration number: "<< count << std::endl;
229
230//adapt initial sublimation to sublimation in next step (take the average)
231 if (count>1){
232 nodes_Subl_ini.grid3D = nodes_Subl_ini.grid3D+nodes_Subl.grid3D;
233nodes_Subl_ini.grid3D = nodes_Subl_ini.grid3D/2.;
234 }
235
236//save ini concentration to compare to next step
237values_elements_to_nodes(nodes_tmp_c,c);
238
239//Combine initial sublimation and initial humidity field and solve new steady state for humidity
240if (Q_FB==1){ //switch humidity feedback
241initializeSystem( colA, rowA,sA,sB,rhs,f,Psi,q,q00, HUM);
242Grid3DObject tmpSource_q;
243tmpSource_q.set(nx,ny,nz,ta.cellsize,ta.llcorner);
244 for (size_t ii=0; ii<(nx*ny*nz); ii++) {
245tmpSource_q(ii) = -1.*nodes_Subl_ini(ii) / (Atmosphere::waterVaporDensity(nodes_Tair(ii), Atmosphere::vaporSaturationPressure(nodes_Tair(ii)))+Atmosphere::stdDryAirDensity(nodes_z(ii), nodes_Tair(ii))*factorf);
246}
247values_nodes_to_elements(tmpSource_q,f);
248prepareSolve();
249SolveEquation(0, maxTimeStep, HUM);
250
251//recalculate relative humidity field
252values_elements_to_nodes(nodes_q, q);
253for (size_t ii=0; ii<(nx*ny*nz); ii++) {
254nodes_RH(ii) = RH_from_q(nodes_Tair(ii),nodes_q(ii), nodes_z(ii));
255if (nodes_RH(ii)>=1.5){
256nodes_q(ii) = Atmosphere::relToSpecHumidity(nodes_z(ii), nodes_Tair(ii),nodes_RH(ii));
257}
258}
259double tmpmin=nodes_q.grid3D.getMin();
260if (tmpmin<0){
261cout<<"WARNING, at least one point with negative specific humidity! (minimum q="<<tmpmin<<")"<<endl;
262}
263}
264
265if (C_FB== 1){
266//recalculate concentration field
267initializeSystem( colA, rowA,sA,sB, rhs, f, Psi, c,c00, CON);
268Grid3DObject tmpSink_c;
269tmpSink_c.set(nx,ny,nz,ta.cellsize,ta.llcorner);
270tmpSink_c.grid3D = nodes_Subl_ini.grid3D*factorf;
271values_nodes_to_elements(tmpSink_c,f);
272prepareSolve();
273timeStep=0;
274SolveEquation(timeStep, maxTimeStep, CON);
275values_elements_to_nodes(nodes_c,c);
276const double tmpmin = nodes_c.grid3D.getMin();
277if (tmpmin<0){
278std::cout<<"!!!WARNING!!! at least one point with negative snow concentration! (minimum c="<<tmpmin*1e3<<" g/m3)"<<std::endl;
279}
280}
281
282//temperature feedback
283if (T_FB==1){
284initializeSystem( colA, rowA,sA,sB,rhs,f,Psi,T,T00, TEM);
285Grid3DObject tmpSink_T;
286tmpSink_T.set(nx,ny,nz,ta.cellsize,ta.llcorner);
287for (size_t ii=0; ii<(nx*ny*nz); ii++) {
288tmpSink_T(ii) = (nodes_Subl_ini(ii)*factorf) * Constants::lh_sublimation * 1/(Constants::specific_heat_air*(1+0.84*nodes_q(ii))*Atmosphere::stdDryAirDensity(nodes_z(ii), nodes_Tair(ii)));
289}
290values_nodes_to_elements(tmpSink_T,f);
291prepareSolve();
292timeStep=0;
293SolveEquation(timeStep, maxTimeStep, TEM);
294
295//get temperature from potential temperature and calculate relative humidity
296Grid3DObject tmp_potT;
297tmp_potT.set(nx,ny,nz,ta.cellsize,ta.llcorner);
298values_elements_to_nodes(tmp_potT,T);
299for (size_t ii=0; ii<(nx*ny*nz); ii++) {
300nodes_Tair(ii) = tmp_potT(ii)-(Cst::gravity/Constants::specific_heat_air)*nodes_z(ii);
301nodes_RH(ii) = RH_from_q(nodes_Tair(ii),nodes_q(ii), nodes_z(ii));
302if (nodes_RH(ii)>=1.5){
303nodes_q(ii)=Atmosphere::relToSpecHumidity(nodes_z(ii), nodes_Tair(ii),nodes_RH(ii));
304}
305}
306const double tmpmin = nodes_Tair.grid3D.getMin();
307if (tmpmin<0){
308std::cout<<"WARNING, at least one point with negative temperature! (minimum T="<<tmpmin<<")"<<std::endl;
309}
310}
311
312Sublimation();
313
314//find maximum difference and calculate mean final sublimation
315Grid3DObject temp_dif_subl;
316Grid3DObject temp_dif_c;
317temp_dif_subl.set(nx,ny,nz,ta.cellsize,ta.llcorner);
318temp_dif_c.set(nx,ny,nz,ta.cellsize,ta.llcorner);
319
320temp_dif_subl.grid3D=(nodes_Subl.grid3D-nodes_Subl_ini.grid3D);
321temp_dif_subl.grid3D.abs();
322
323temp_dif_c.grid3D=nodes_c.grid3D-nodes_tmp_c.grid3D;
324temp_dif_c.grid3D.abs();
325
326maxdif_s=temp_dif_subl.grid3D.getMax();
327maxdif_c=temp_dif_c.grid3D.getMax();
328mean_endsubl=nodes_Subl.grid3D.getMean();
329
330//test whether changes in c and sublimation are too large
331if (C_FB==1){
332 if (maxdif_s>1.2e-6 && maxdif_c>1.1e-6){
333test=0; //go on with iteration
334 } else {
335test=1; //found steady state, stop iteration
336 }
337} else {//test without concentration as this will not change when feedback is switched off
338 if (maxdif_s>1.e-6){
339test=0;
340 }else{
341test=1;
342 }
343}
344std::cout <<"Maximum difference in sublimation= " << maxdif_s*1e6 << "e-6, maximum difference in c= " << maxdif_c*1e3 << "e-3"<<std::endl;
345}//close while
346
347std::cout<<"Mean final sublimation=" <<mean_endsubl<<", fraction of initial:"<<mean_endsubl/(mean_startsubl+1e-20)<<std::endl;
348
349if (test==1){
350 std::cout <<"Found steady state sublimation after "<<count<<" iterations."<< std::endl;
351}else{
352 std::cout<<"!!!WARNING!!! No steady state found yet (maximum difference in sublimation = "<<maxdif_s<<"), max number of iterations reached, using sublimation field anyway"<<std::endl;
353}
354
355}
356
357/**
358 * @brief prepareSolve
359 * Some preparations to solve a 3D-field.
360 */
361void SnowDriftA3D::prepareSolve(){
362
363 prepareSparseMatrix( colA,rowA,adjA);
364 resetArray( sA );
365 resetArray( sB );
366 resetArray( rhs );
367 assembleSystem( colA, rowA,sA,sB,Psi,f,dt_diff);
368 //-------------------------------------------------------------------------------
369 // Apply BC-----------NOT implemented!!!!!!!!!!!!!!!!!!!!!!!
370 //-------------------------------------------------------------------------------
371 //applyBoundaryValues(c00,Psi);
372}
373
374

Archive Download this file

Revision: HEAD