Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/snowdrift/SnowDriftFENumerics.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#include <ctime>
23
24#include <meteoio/MeteoIO.h>
25
26#include <alpine3d/snowdrift/SnowDrift.h>
27
28using namespace std;
29
30//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31//
32// Function : setQuadraturePoints
33// Authors : Henning Loewe
34//--------------------------------------------------------------------------------
35// Description : sets the SnowDrift variable qValues which holds the
36// quadrature points for a one point gaussian quadrature
37//
38// IP_8-----IP_7
39// /| /| /\ z
40// / | / | |
41// / | / | |
42// IP_5-----IP_6 | |------> y
43// | | | | /
44// | IP_4--|--IP_3 /
45// | / | / /
46// | / | / V x
47// |/ |/
48// IP_1------IP_2
49//--------------------------------------------------------------------------------
50//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51void SnowDriftA3D::setQuadraturePoints( void )
52{
53 double IP[3][8]={ { 1, 1,-1,-1, 1, 1,-1,-1},
54 {-1, 1, 1,-1,-1, 1, 1,-1},
55 {-1,-1,-1,-1, 1, 1, 1, 1} };
56
57 double qValue = 0.577350269189626;
58 for ( int i = 0; i < 3; i++){
59 for ( int j = 0; j < 8; j++){
60 qPoint(i,j) = qValue * IP[i][j];
61}
62 }
63}
64
65/**
66 * @brief Compute Jacobian
67 * This function calculates the Jacobian at the given
68 * point P and and stores the value of its determinant at point P in
69 * the global variable DETERMINANTJ[k] the input value is the
70 * respective element, the position where to store the
71 * DETERMINANTJ[k], and the function is of type void, the Jacobian
72 * matrix J and the integration point P and the DETERMINANTJ array
73 * @param DETERMINANTJ Determinant of Jacobian matrix
74 * @param J Jacobian matrix
75 * @param element element number
76 * @param P point where Jacobian has to be calculated
77 * @param k for storing Jacobian at DETERMINANTJ(k)
78 */
79void SnowDriftA3D::Jacobian(double *DETERMINANTJ,
80 double J[][3],
81 const int element,
82 const double *P,
83 const int k,
84 const int ix,
85 const int iy,
86 const int iz)
87{
88 double G[3][8]; //dummy array
89
90 // the partial derivatives of phi evaluated at the integration point P
91 TTfun(G,P);
92 // set Jacobian J to zero
93 for (int i = 0; i < 3; i++) {
94 for ( int j = 0; j < 3; j++) {
95 J[i][j] = 0;
96}
97 }
98
99 int ii = 0, jj = 0, kk = 0;
100
101 for (int j=0; j<3; j++){
102int count=0;
103while (count<8){
104 if (count==0 || count==3 || count==4 || count==7){
105ii=ix;
106 } else{
107ii=ix+1;
108 }
109 if (count==0 || count==1 || count==4 || count==5){
110jj=iy;
111 } else{
112jj=iy+1;
113 }
114 if (count<=3){
115kk=iz;
116 } else{
117kk=iz+1;
118 }
119
120 J[0][j] = J[0][j] + nodes_x.grid3D(ii,jj,kk) * G[j][count];
121 J[1][j] = J[1][j] + nodes_y.grid3D(ii,jj,kk) * G[j][count];
122 J[2][j] = J[2][j] + nodes_z.grid3D(ii,jj,kk) * G[j][count];
123 count++;
124}
125 }
126
127
128 // now calculate the determinant
129 DETERMINANTJ[k] = J[0][0] * (J[1][1]*J[2][2]-J[2][1]*J[1][2])
130 -J[0][1]*(J[1][0]*J[2][2]-J[2][0]*J[1][2])
131 +J[0][2]*(J[1][0]*J[2][1]-J[2][0]*J[1][1]);
132
133 if (DETERMINANTJ[k]<0){
134 cout<<"the determinant of the element "<<element<<" is negative, watch out"<<endl;
135 }
136
137}
138
139/**
140 * @brief J0 function
141 * computes the matrix of cofactors of the input matrix J and stores it in J0
142 * @param J0 matrix of cofactors (transpose of the adjoint of J)
143 * @param J Jacobian matrix
144 */
145void SnowDriftA3D::J0fun( double J0[3][3], const double J[3][3])
146{
147 J0[0][0] = J[1][1]*J[2][2]-J[2][1]*J[1][2];
148 J0[0][1] = -(J[1][0]*J[2][2]-J[2][0]*J[1][2]);
149 J0[0][2] = J[1][0]*J[2][1]-J[1][1]*J[2][0];
150
151 J0[1][0] = -(J[0][1]*J[2][2]-J[2][1]*J[0][2]);
152 J0[1][1] = J[0][0]*J[2][2]-J[2][0]*J[0][2];
153 J0[1][2] = -(J[0][0]*J[2][1]-J[2][0]*J[0][1]);
154
155 J0[2][0] = J[0][1]*J[1][2]-J[1][1]*J[0][2];
156 J0[2][1] = -(J[0][0]*J[1][2]-J[1][0]*J[0][2]);
157 J0[2][2] = J[0][0]*J[1][1]-J[1][0]*J[0][1];
158}
159
160/**
161 * @brief GQIntB
162 * returns the integral of the product of basis
163 * functions i and j over a hexahedral element using gaussian
164 * quadrature, DETERMINANTK[0],...DETERMINANTK[7] contains the
165 * determinant of the jacobian of the mapping onto the unit cube
166 * evaluated at the corners
167 * @param DETERMINANTJ determinant of Jacobian
168 * @param i nr shape function
169 * @param j nr shape function
170 */
171double SnowDriftA3D::GQIntB( double *DETERMINANTJ,const int i, const int j)
172{
173 double q_loc=0;
174 double PHI[8];
175 int k;
176 double cad[3];
177
178 for ( k = 0; k < 8; k++) {
179 cad[0]=qPoint(0,k);
180 cad[1]=qPoint(1,k);
181 cad[2]=qPoint(2,k);
182
183 phi(PHI,cad);
184
185 q_loc += PHI[i] * PHI[j] * fabs(DETERMINANTJ[k]); // i and j have already been adjustet!
186 }
187 return q_loc;
188}
189
190/**
191 * @brief GQIntC
192 * returns the integral of the product of the gradients
193 * of basis functions i and j over a hexahedral element using gaussian
194 * quadrature, DETERMINANTK[0],...DETERMINANTK[7] contains the
195 * determinant of the jacobian of the mapping onto the unit cube
196 * evaluated at the corners
197 * @param DETERMINANTJ determinant of Jacobian
198 * @param J0M
199 * @param i pointer/index
200 * @param j pointer/index
201 * @param b drift vector
202 * @param K Diffusion tensor
203 */
204double SnowDriftA3D::GQIntC(double * DETERMINANTJ,
205 const double J0M[3][3][8],
206 const int i, const int j,
207 const double b[3],
208 const double K[3][3])
209{
210 double q_loc=0;
211 double PHI[8], TT[3][8], Hi[3], Hj[3];
212 double A[3], B[3], C[3]; //pure dummy vectors
213 int k,m;
214 double cad[3];
215
216
217 for ( k = 0;k<8;k++){
218
219 cad[0]=qPoint(0,k);//the integration point
220 cad[1]=qPoint(1,k);
221 cad[2]=qPoint(2,k);
222
223 phi(PHI,cad);
224 TTfun(TT,cad);
225
226 Hi[0]=TT[0][i];//dummy vector containing the gradient of the ith shape function
227 Hi[1]=TT[1][i];
228 Hi[2]=TT[2][i];
229
230 Hj[0]=TT[0][j];// dummy vector containing the gradient of the jth shape function
231 Hj[1]=TT[1][j];
232 Hj[2]=TT[2][j];
233
234 /* the first part is due to diffusion */
235
236 /* first calculate the intermediate vectors*/
237
238 for (m=0;m<3;m++){ /* set to zero*/
239 A[m]=0;
240 B[m]=0;
241 C[m]=0;
242 }
243
244 for (m=0;m<3;m++){
245 A[0]=A[0]+J0M[0][m][k]*Hj[m]; /* this is J0M*Hj */
246 A[1]=A[1]+J0M[1][m][k]*Hj[m];
247 A[2]=A[2]+J0M[2][m][k]*Hj[m];
248 }
249
250 for (m=0;m<3;m++){
251 B[0]=B[0]+K[0][m]*A[m]; /* this is K*A */
252 B[1]=B[1]+K[1][m]*A[m];
253 B[2]=B[2]+K[2][m]*A[m];
254
255 }
256
257 for (m=0;m<3;m++){
258 C[0]=C[0]+J0M[0][m][k]*Hi[m];/* this is J0M*Hi */
259 C[1]=C[1]+J0M[1][m][k]*Hi[m];
260 C[2]=C[2]+J0M[2][m][k]*Hi[m];
261 }
262
263
264 for (m=0;m<3;m++) {
265 q_loc=q_loc+B[m]*C[m]*(1/fabs(DETERMINANTJ[k])); /* add up the diffusion part */
266}
267
268 /* now we add the part that is due to convection*/
269
270 // A is the matrix product J0M(k)*Hj and has been calculated previously, no need to repeat this calculation
271
272 if (fabs(DETERMINANTJ[k])<1e-10)
273 printf("\n the determinant of J at point k is zero or almost zero, watch out\n");
274
275 for (m=0;m<3;m++){
276 q_loc=q_loc+b[m]*A[m]*PHI[i]*(DETERMINANTJ[k]/fabs(DETERMINANTJ[k]));
277}
278 }//end loop in k
279
280 return q_loc;
281
282}
283
284/**
285 * @brief GQIntApdx
286 * returns the integral of the product of a basis
287 * function i and the gradient of basis function j over a hexahedral
288 * element using gaussian quadrature,
289 * DETERMINANTK[0],...DETERMINANTK[7] contains the determinant of the
290 * jacobian of the mapping onto the unit cube evaluated at the corners
291 * @param DETERMINANTJ determinant of Jacobian
292 * @param J0M
293 * @param i nr basis function
294 * @param j nr basis function
295 * @param b drift vector
296 * @param deltak parameter SUPG method
297 */
298double SnowDriftA3D::GQIntApdx(double *DETERMINANTJ,
299 const double J0M[3][3][8],
300 const int i, const int j,
301 double b[],
302 const double deltak)
303{
304 double q_loc=0;
305 double A[3];
306 double PHI[8], TT[3][8], Hi[3];
307 double cad[3];
308 int m,k;
309
310 for (k=0;k<8;k++){
311
312 cad[0]=qPoint(0,k);
313 cad[1]=qPoint(1,k);
314 cad[2]=qPoint(2,k);
315
316 phi(PHI,cad);
317 TTfun(TT,cad);
318
319 Hi[0]=TT[0][i];//dummy vector containing the gradient of the ith shape function
320 Hi[1]=TT[1][i];
321 Hi[2]=TT[2][i];
322
323 for (m=0;m<3;m++){
324 A[m]=0;
325}
326
327 for (m=0;m<3;m++){
328 A[0]+=J0M[0][m][k]*Hi[m];
329 A[1]+=J0M[1][m][k]*Hi[m];
330 A[2]+=J0M[2][m][k]*Hi[m];
331 }
332
333
334 for (m=0;m<3;m++){
335 q_loc+=PHI[j]*b[m]*A[m]*(DETERMINANTJ[k]/fabs(DETERMINANTJ[k]));
336 }
337}
338 q_loc=q_loc*deltak;
339 return q_loc;
340}
341
342/**
343 * @brief GQIntAdxdx
344 * ...
345 * @param DETERMINANTJ determinant of Jacobian
346 * @param J0M
347 * @param i pointer
348 * @param j pointer
349 * @param *b drift vector
350 * @param deltak parameter SUPG method
351 */
352double SnowDriftA3D::GQIntAdxdx( double *DETERMINANTJ,
353 double J0M[3][3][8],
354 const int i, const int j,
355 double *b,
356 const double deltak)
357{
358 double q_loc=0; double z_loc,y_loc;
359 double TT[3][8], Hi[3], Hj[3];
360 double A[3],C[3];
361 double cad[3];
362
363 for (int k = 0; k < 8; k++){
364
365 cad[0] = qPoint(0,k); //the kth integration point
366 cad[1] = qPoint(1,k);
367 cad[2] = qPoint(2,k);
368
369 TTfun(TT,cad);
370
371 Hi[0] = TT[0][i];//dummy vector containing the gradient of the ith shape function
372 Hi[1] = TT[1][i];
373 Hi[2] = TT[2][i];
374
375 Hj[0] = TT[0][j];// dummy vector containing the gradient of the jth shape function
376 Hj[1] = TT[1][j];
377 Hj[2] = TT[2][j];
378
379 /* the first part is due to diffusion */
380
381 /* first calculate the intermediate vectors*/
382 for (int m = 0; m < 3; m++){
383 A[m] = 0;
384 C[m] = 0;
385}
386 z_loc = 0;
387 y_loc = 0;
388
389 for (int m = 0; m < 3; m++){
390 A[0] = A[0] + J0M[0][m][k] * Hj[m]; /* this is J0M*Hj */
391 A[1] = A[1] + J0M[1][m][k] * Hj[m];
392 A[2] = A[2] + J0M[2][m][k] * Hj[m];
393}
394
395 for (int m = 0;m < 3; m++){
396 C[0]=C[0]+J0M[0][m][k]*Hi[m];/* this is J0M*Hi */
397 C[1]=C[1]+J0M[1][m][k]*Hi[m];
398 C[2]=C[2]+J0M[2][m][k]*Hi[m];
399}
400
401 for (int m = 0;m < 3; m++){
402 z_loc+=b[m]*A[m];
403 y_loc+=b[m]*C[m];
404}
405
406 q_loc+=z_loc*y_loc*(1/fabs(DETERMINANTJ[k])); /* add up the diffusion part */
407 }
408
409 q_loc = q_loc*deltak;
410 return q_loc;
411
412}
413
414/**
415 * @brief TT function
416 * TTfun is a function calculating the values of the
417 * partial derivative of phi at the point P and putting them in a
418 * (3,8) array
419 * @param TT partial derivative of phi at point P (result of this function)
420 * @param P pointer
421 */
422void SnowDriftA3D::TTfun(double TT[3][8],const double *P)// P is a pointer pointing at an adress containing a double
423{
424 const double p_loc = P[0];
425 const double q_loc = P[1];
426 const double r = P[2];//the three coordiantes of the point P
427
428 TT[0][0] = 1/8.*(1-q_loc)*(1-r);
429 TT[0][1] = 1/8.*(1+q_loc)*(1-r);
430 TT[0][2] = -1/8.*(1+q_loc)*(1-r);
431 TT[0][3] = -1/8.*(1-q_loc)*(1-r);
432 TT[0][4] = 1/8.*(1-q_loc)*(1+r);
433 TT[0][5] = 1/8.*(1+q_loc)*(1+r);
434 TT[0][6] = -1/8.*(1+q_loc)*(1+r);
435 TT[0][7] = -1/8.*(1-q_loc)*(1+r);
436
437 TT[1][0] = -1/8.*(1+p_loc)*(1-r);
438 TT[1][1] = 1/8.*(1+p_loc)*(1-r);
439 TT[1][2] = 1/8.*(1-p_loc)*(1-r);
440 TT[1][3] = -1/8.*(1-p_loc)*(1-r);
441 TT[1][4] = -1/8.*(1+p_loc)*(1+r);
442 TT[1][5] = 1/8.*(1+p_loc)*(1+r);
443 TT[1][6] = 1/8.*(1-p_loc)*(1+r);
444 TT[1][7] = -1/8.*(1-p_loc)*(1+r);
445
446 TT[2][0] = -1/8.*(1+p_loc)*(1-q_loc);
447 TT[2][1] = -1/8.*(1+p_loc)*(1+q_loc);
448 TT[2][2] = -1/8.*(1-p_loc)*(1+q_loc);
449 TT[2][3] = -1/8.*(1-p_loc)*(1-q_loc);
450 TT[2][4] = 1/8.*(1+p_loc)*(1-q_loc);
451 TT[2][5] = 1/8.*(1+p_loc)*(1+q_loc);
452 TT[2][6] = 1/8.*(1-p_loc)*(1+q_loc);
453 TT[2][7] = 1/8.*(1-p_loc)*(1-q_loc);
454
455}
456
457/**
458 * @brief Phi - shape functions
459 * Calculates the value of the 8 shape functions at point P
460 * @param *PHI ...................shape functions
461 * @param *P pointer
462 */
463void SnowDriftA3D::phi(double *PHI,double *P)
464{
465 double p_loc=P[0], q_loc=P[1], r=P[2];
466
467 PHI[0] = 1/8.*(1+p_loc)*(1-q_loc)*(1-r);
468 PHI[1] = 1/8.*(1+p_loc)*(1+q_loc)*(1-r);
469 PHI[2] = 1/8.*(1-p_loc)*(1+q_loc)*(1-r);
470 PHI[3] = 1/8.*(1-p_loc)*(1-q_loc)*(1-r);
471 PHI[4] = 1/8.*(1+p_loc)*(1-q_loc)*(1+r);
472 PHI[5] = 1/8.*(1+p_loc)*(1+q_loc)*(1+r);
473 PHI[6] = 1/8.*(1-p_loc)*(1+q_loc)*(1+r);
474 PHI[7] = 1/8.*(1-p_loc)*(1-q_loc)*(1+r);
475}
476
477/**
478 * @brief matmult
479 * computes a matrix vector product for a sparse matrix
480 * @param res result (matrix vector product)
481 * @param x
482 * @param sm sparse matrix
483 * @param ijm
484 */
485void SnowDriftA3D::matmult(CDoubleArray& res, const CDoubleArray& x_loc, double* sm, int* ijm)
486{
487 int n= ijm[1]-2;
488
489 for (int i=1;i<=n;i++){
490 res[i]=0;
491}
492
493 for (int i = 1; i <= n; i++){
494 res[i]=sm[i]*x_loc[i];
495 for (int k = ijm[i]; k <= (ijm[i+1]-1); k++) {
496 //res[i]=ijm[5];//fake line
497 //res[i]=x[ijm[k]];//fake line
498 res[i] = res[i] + sm[k]*x_loc[ijm[k]]; //real line
499}
500 }
501}
502
503/**
504 * @brief matmult
505 * computes a matrix vector product for a sparse matrix of CSR format
506 * @param y result (matrix vector product)
507 * @param x
508 * @param sA
509 * @param colind column index
510 * @param rowPtr row index
511 */
512void SnowDriftA3D::matmult(CDoubleArray& y_loc,
513const CDoubleArray& x_loc,
514const CDoubleArray& sA_loc,
515const CIntArray& colInd,
516CIntArray& rowPtr )
517{
518const size_t dim = rowPtr.getNx() - 1;
519
520for (size_t i = 0; i < dim; i++) {
521y_loc[i] = 0;
522for (int j = rowPtr[i]; j < rowPtr[i+1] ; j++) {
523y_loc[i] += sA_loc[ j ] * x_loc[ colInd[j] ];
524}
525}
526}
527
528/**
529 * @brief transmult
530 * this is the transmult function that multiplies the vector x by the transpose of the matrix M in its
531 * sparse form sm, ijm
532 * @param res result
533 * @param x
534 * @param sm
535 * @param ijm index
536 */
537void SnowDriftA3D::transmult(CDoubleArray& res, const CDoubleArray& x_loc, double* sm, int* ijm)
538{
539 int n=ijm[1]-2;
540
541 int j;
542 for (int i = 1; i <= n; i++) {
543 res[i] = sm[i]*x_loc[i];
544 }
545
546 for (int i = 1; i<= n;i++) {
547 for (int k = ijm[i]; k <= (ijm[i+1]-1); k++) {
548 j = ijm[k];
549 res[j] = res[j] + sm[k] * x_loc[i];
550}
551 }
552}
553
554/**
555 * @brief bicgStab iterative equation solver
556 * iterative equation solver
557 * Tests : Tested by the followin procedure: given a sparse matrix A
558 * (CRS-format) characterized by colA and rowA and an arbitrary, or
559 * rather: a few nontrivial examples of a vector x. For each x
560 * matmult(y,x,sA,rowA,colA) and bicgStab(result,y,sA,colA,rowA,...)
561 * have been computed and then verified that result=x
562 * @param res result
563 * @param rhs
564 * @param sA
565 * @param colA
566 * @param rowA
567 * @param nmax max number of iterations
568 * @param tol tolerance
569 */
570void SnowDriftA3D::bicgStab(CDoubleArray& result,
571 CDoubleArray& rhs_loc,
572 const CDoubleArray& sA_loc,
573 const CIntArray& colA_loc,
574 CIntArray& rowA_loc,
575 const int nmax,
576 const double tol,
577 double& testres)
578{
579
580 //dimension of the system
581 const size_t n = rowA_loc.getNx() - 1;
582
583 CDoubleArray res1;
584 res1.resize( n);
585
586 CDoubleArray r_0;
587 r_0.resize( n );
588
589 CDoubleArray r;
590 r.resize( n );
591
592 CDoubleArray p_loc;
593 p_loc.resize( n );
594
595 CDoubleArray v;
596 v.resize( n );
597
598 CDoubleArray aux1;
599 aux1.resize( n );
600
601 CDoubleArray auxhat;
602 auxhat.resize( n );
603
604 CDoubleArray phat;
605 phat.resize( n );
606
607 CDoubleArray aux2;
608 aux2.resize( n );
609
610 double rho_old = 1;
611 double rho_new = 1;
612
613 double omega = 1;
614 double alpha = 1;
615 double beta = 0;
616
617 double residual,norm1,norm2;
618 double res4,res5;
619 double tmp_res=1;
620 int iterations=0;
621
622 CDoubleArray res_ok_state;
623
624
625 srand(static_cast<unsigned>(time(NULL)));
626 //intitialization
627 for (size_t i=0;i<n;i++) {
628 result[i]=0;//(rand()%RAND_MAX)/(1.0*RAND_MAX);
629 v[i] = 0;
630 p_loc[i] = 0;
631 }
632
633 matmult(res1,result,sA_loc,colA_loc,rowA_loc); // multiply Bx and store it into the dummy res1
634
635 norm1 = 0;
636 norm2 = 0;
637 for ( size_t i = 0; i < n; i++ ) {
638 r_0[i] = rhs_loc[i]-res1[i];
639//cout<<"r_0="<<r_0[i]<<" rhs= "<<rhs[i]<<" res1="<<res1[i]<<"at i"<<i<<endl;
640
641 r[i] = rhs_loc[i]-res1[i];
642 norm1 += rhs_loc[i] * rhs_loc[i];
643 norm2 += r_0[i] * r_0[i];
644 }
645
646 if ( rho_new == 0 ){
647 printf("Iteration fails\n"); //in this case, can't take x=0 as initial guess
648}
649
650 int k = 0;
651 double mark = 0; //as soon as mark==1 you can stop the iteration, good approximatin is attained
652
653 residual = sqrt(norm2);
654
655 if ( residual < tol) {//stopping criterion
656 mark=1;
657 printf("-------> Lucky failure within %d steps with residual: %f\n", k, residual);
658
659 }
660
661 //main loop
662 while ( (k<=nmax) && (mark==0) ) {
663
664 rho_new = 0;
665 for (size_t i=0;i<n;i++){
666 rho_new += r_0[i]*r[i];
667}
668
669 beta = rho_new / rho_old * alpha / omega;
670
671 for (size_t i=0;i<n;i++) {
672 p_loc[i] *= beta;
673 p_loc[i] += (r[i] - beta * omega *v[i] );
674
675 //precond
676 phat[i] = p_loc[i]/precond[i];
677}
678
679 matmult(v,phat,sA_loc,colA_loc,rowA_loc);// put B*p into v
680
681 res4 = 0;
682 for ( size_t i = 0; i < n ; i++ ){
683 res4 += v[i] * r_0[i];
684}
685
686 //HERE test print if v or res4 don't make sense anymore
687 if (!(fabs(v[n-2])<1e20)){
688printf("-------> LH: v too large?\n");
689}
690 if (!(fabs(res4)<1e40)){
691printf("-------> LH: res4 too large? res4=%f \n", res4);
692}
693 alpha = rho_new / res4;
694
695 for ( size_t i = 0; i < n; i++ ) {
696 aux1[i] = r[i] - alpha * v[i];
697 //precond
698 auxhat[i] = aux1[i]/precond[i];
699 }
700
701 matmult(aux2,auxhat,sA_loc,colA_loc,rowA_loc);
702
703 res4 = 0;
704 res5 = 0;
705 for ( size_t i = 0; i < n; i++ ) {
706 res4 += ( aux2[i] * aux1[i] );
707 res5 += ( aux2[i] * aux2[i] );
708 }
709
710 omega = res4 / res5;
711
712 for ( size_t i = 0; i < n; i++ ) {
713 result[i] += ( alpha * phat[i] + omega * auxhat[i] );
714 r[i] = aux1[i] - omega * aux2[i];
715 }
716
717 matmult(res1,result,sA_loc,colA_loc,rowA_loc);
718
719 norm1=0;
720 norm2=0;
721 for ( size_t i = 0; i < n; i++ ) {
722 norm1 += ( rhs_loc[i] * rhs_loc[i] );
723norm2 +=((rhs_loc[i]-res1[i])*(rhs_loc[i]-res1[i])) ; //norm2 += pow( rhs[i] - res1[i], 2 );
724 }
725
726 residual = sqrt(norm2) / sqrt(norm1);
727
728 if ( residual <= tol){//stopping criteria!
729 mark=1;
730 printf("-------> Convergence within %d steps ", k);
731 printf("with residual: %f \n", residual);
732 }
733
734 if ( residual >= 1.e4 && tmp_res>=1.) {//stopping criteria!
735 mark=1;
736 printf("-------> Hopeless after %d steps ", k);
737 printf("with residual: %f, try again! \n ", residual);
738 }
739
740 if ( residual < tmp_res){
741 //copy this state temporarily
742 res_ok_state=result;
743 tmp_res=residual;
744 iterations=k;
745 }
746
747 rho_old = rho_new;
748// //LH-CHANGES BEGIN
749// // 1) comment out the following line
750// cout<<" this is k: "<<k<< " and the residual: "<< residual<<endl;
751// // 2) redirect output to logfile
752// FILE *logfile;
753// string fname=string("../results/numerics-tests/residual.dat");
754// logfile = fopen(fname.c_str(), "a");
755// fprintf(logfile, "%d %.15lf\n", k, residual);
756// fclose (logfile);
757// //LH-CHANGES END
758 k++;
759 }
760 if (mark==0) {
761printf("-------> No convergence within maximal number of iterations (residual = %f)\n",residual);
762if (residual > 1e7*tol && tmp_res<1.){
763 printf("Use a previous step with residual = %f, #iterations= %d)\n", tmp_res, iterations);
764 //copy previous result and residual
765 result=res_ok_state;
766 residual=tmp_res;
767 }
768 }
769 testres=residual;
770
771}//end of function

Archive Download this file

Revision: HEAD