[44] | 1 | /*************************************************************************** |
---|
| 2 | module classe solsor_gc1_x_r.h - description |
---|
| 3 | ***************************************************************************/ |
---|
| 4 | // Julien Brajard |
---|
| 5 | // locean-ipsl.upmc, Paris, April 25, 2011 |
---|
| 6 | //=========================================================================== |
---|
| 7 | // methode forward |
---|
| 8 | forward (YREAL x1, YREAL x2, YREAL x3, YREAL x4, YREAL x5, YREAL x6, YREAL x7, YREAL x8) |
---|
| 9 | { |
---|
| 10 | /* 1 from gcx_dynspg_flt 1 i j t |
---|
| 11 | 2 from gcb_dynspg_flt 1 i j t |
---|
| 12 | 3 from solsor_gc2_x_r 1 i j-1 k-1 t |
---|
| 13 | 4 from solsor_gc2_x_r 1 i-1 j k-1 t |
---|
| 14 | 5 from solsor_gc2_x_r 1 i+1 j k-1 t |
---|
| 15 | 6 from solsor_gc2_x_r 1 i j+1 k-1 t |
---|
| 16 | 7 from solsor_gc2_x_r 1 i j k-1 t |
---|
| 17 | 8 from solsor_flag 2 NX NY k-1 t |
---|
| 18 | } |
---|
| 19 | output 1 : gcx |
---|
| 20 | oupput 2 : gcr (pour l'arrêt des itérations |
---|
| 21 | |
---|
| 22 | */ |
---|
[45] | 23 | if(Yt==TU) |
---|
| 24 | { |
---|
| 25 | YS1=0; |
---|
| 26 | YS2=0; |
---|
| 27 | } |
---|
| 28 | else { |
---|
| 29 | int t; |
---|
| 30 | if(Yt==TU+1 && neuler==0) |
---|
| 31 | t=0; |
---|
| 32 | else |
---|
| 33 | t=1; |
---|
| 34 | |
---|
| 35 | if (x8==1) |
---|
| 36 | { |
---|
| 37 | YS1=x7; |
---|
| 38 | YS2=0; |
---|
| 39 | } |
---|
| 40 | else |
---|
| 41 | { |
---|
| 42 | if (Yk==0) { |
---|
| 43 | if ( ( (Yi==0) | (Yj==0) ) | ( (Yi==NX-1) | (Yj==NY-1) )) { |
---|
| 44 | YS1=0; |
---|
| 45 | YS2=0; |
---|
| 46 | } |
---|
| 47 | else { |
---|
| 48 | YS1=x1; |
---|
| 49 | YS2=0; |
---|
| 50 | } |
---|
| 51 | } //Yk=0; |
---|
| 52 | else { |
---|
| 53 | YREAL ztmp,zres; |
---|
| 54 | int ishift=Yj%2; |
---|
| 55 | if ((Yi+1)%2==ishift) { //! Guess black update |
---|
| 56 | ztmp = x2 |
---|
| 57 | - gcp(Yi,Yj,0,t) * x3 |
---|
| 58 | - gcp(Yi,Yj,1,t) * x4 |
---|
| 59 | - gcp(Yi,Yj,2,t) * x5 |
---|
| 60 | - gcp(Yi,Yj,3,t) * x6; |
---|
| 61 | // ! Estimate of the residual |
---|
| 62 | zres = ztmp - x7; |
---|
| 63 | YS2 = zres * gcdmat(Yi,Yj,t) * zres; |
---|
| 64 | // ! Guess update |
---|
| 65 | YS1 = sor * ztmp + (1.-sor) * x7; |
---|
| 66 | } |
---|
| 67 | else |
---|
| 68 | { |
---|
| 69 | YS2=0; |
---|
| 70 | YS1=x7; |
---|
| 71 | } |
---|
| 72 | |
---|
| 73 | |
---|
| 74 | } // Yk>0 |
---|
| 75 | } //if x8==0 |
---|
[49] | 76 | //if (Yt==TU+1) |
---|
| 77 | // printf("%d\t%d\t%d\t%e\t%e (1)\n",Yk,Yj,Yi,YS1,YS2); |
---|
[45] | 78 | } //Yt>TU |
---|
[44] | 79 | } |
---|
| 80 | |
---|
| 81 | //=========================================================================== |
---|
| 82 | // methode backward |
---|
| 83 | |
---|
[45] | 84 | backward (YREAL x1, YREAL x2, YREAL x3, YREAL x4, YREAL x5, YREAL x6, YREAL x7, YREAL x8) { |
---|
[49] | 85 | /* 1 from gcx_dynspg_flt 1 i j t |
---|
| 86 | 2 from gcb_dynspg_flt 1 i j t |
---|
| 87 | 3 from solsor_gc2_x_r 1 i j-1 k-1 t |
---|
| 88 | 4 from solsor_gc2_x_r 1 i-1 j k-1 t |
---|
| 89 | 5 from solsor_gc2_x_r 1 i+1 j k-1 t |
---|
| 90 | 6 from solsor_gc2_x_r 1 i j+1 k-1 t |
---|
| 91 | 7 from solsor_gc2_x_r 1 i j k-1 t |
---|
| 92 | 8 from solsor_flag 2 NX NY k-1 t |
---|
| 93 | } |
---|
| 94 | output 1 : gcx |
---|
| 95 | oupput 2 : gcr (pour l'arrêt des itérations |
---|
| 96 | |
---|
| 97 | */ |
---|
| 98 | if(Yt==TU) |
---|
| 99 | { |
---|
| 100 | YJ1I1=0; |
---|
| 101 | } |
---|
| 102 | else { |
---|
| 103 | int t; |
---|
| 104 | if(Yt==TU+1 && neuler==0) |
---|
| 105 | t=0; |
---|
| 106 | else |
---|
| 107 | t=1; |
---|
| 108 | |
---|
| 109 | if (x8==1) |
---|
| 110 | { |
---|
| 111 | YJ1I7=1; |
---|
| 112 | // YS1=x7; |
---|
| 113 | } |
---|
| 114 | else |
---|
| 115 | { |
---|
| 116 | if (Yk==0) { |
---|
| 117 | if ( ( (Yi==0) | (Yj==0) ) | ( (Yi==NX-1) | (Yj==NY-1) )) { |
---|
| 118 | |
---|
| 119 | //YS1=0; |
---|
| 120 | |
---|
| 121 | } |
---|
| 122 | else { |
---|
| 123 | YJ1I1=1; |
---|
| 124 | // YS1=x1; |
---|
| 125 | } |
---|
| 126 | } //Yk=0; |
---|
| 127 | else { |
---|
| 128 | // YREAL ztmp,zres; |
---|
| 129 | int ishift=Yj%2; |
---|
| 130 | if ((Yi+1)%2==ishift) { //! Guess black update |
---|
| 131 | YREAL dztmpdx2=1; |
---|
| 132 | YREAL dztmpdx3=-gcp(Yi,Yj,0,t); |
---|
| 133 | YREAL dztmpdx4=-gcp(Yi,Yj,1,t); |
---|
| 134 | YREAL dztmpdx5=-gcp(Yi,Yj,2,t); |
---|
| 135 | YREAL dztmpdx6=-gcp(Yi,Yj,3,t); |
---|
| 136 | /*ztmp = x2 |
---|
| 137 | - gcp(Yi,Yj,0,t) * x3 |
---|
| 138 | - gcp(Yi,Yj,1,t) * x4 |
---|
| 139 | - gcp(Yi,Yj,2,t) * x5 |
---|
| 140 | - gcp(Yi,Yj,3,t) * x6;*/ |
---|
| 141 | |
---|
| 142 | YJ1I2=sor*dztmpdx2; |
---|
| 143 | YJ1I3=sor*dztmpdx3; |
---|
| 144 | YJ1I4=sor*dztmpdx4; |
---|
| 145 | YJ1I5=sor*dztmpdx5; |
---|
| 146 | YJ1I6=sor*dztmpdx6; |
---|
| 147 | YJ1I7=(1.-sor); |
---|
| 148 | //YS1 = sor * ztmp + (1.-sor) * x7; |
---|
| 149 | } |
---|
| 150 | else |
---|
| 151 | { |
---|
| 152 | YJ1I7=1; |
---|
| 153 | // YS1=x7; |
---|
| 154 | } |
---|
| 155 | |
---|
| 156 | |
---|
| 157 | } // Yk>0 |
---|
| 158 | } //if x8==0 |
---|
| 159 | //if (Yt==TU+1) |
---|
| 160 | // printf("%d\t%d\t%d\t%e\t%e (1)\n",Yk,Yj,Yi,YS1,YS2); |
---|
| 161 | } //Yt>TU |
---|
| 162 | |
---|
| 163 | /*if (Yi>22 & Yi < 26 & Yj>18 & Yj <21 & Yk==241 & Yt==9) { |
---|
| 164 | |
---|
| 165 | for (int i=0;i< 9;i++) |
---|
| 166 | { |
---|
| 167 | |
---|
| 168 | for (int j=0;j<2;j++) { |
---|
| 169 | // if (Yjac[j][i]>1 & YG1Y_solsor_gx2_x_r>1e-10 ) { |
---|
| 170 | printf("gx1 (%d,%d,%d,%d) YG1=%e,YG2=%e(%e)",Yi,Yj,Yk,Yt,YG1Y_solsor_gx2_x_r,YG2Y_solsor_gx2_x_r,(&YG1Y_solsor_gx2_x_r)[1]); |
---|
| 171 | printf(" YJ%dI%d=%e\n",j+1,i+1,Yjac[j][i]); |
---|
| 172 | } |
---|
| 173 | } |
---|
| 174 | }*/ |
---|
[44] | 175 | } |
---|