source: trunk/src/solsor_gx1_x_r.h @ 54

Last change on this file since 54 was 49, checked in by jbrlod, 13 years ago

reorganisation de order

  • Property svn:eol-style set to native
File size: 4.0 KB
Line 
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
8forward (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}
19output 1 : gcx
20oupput 2 : gcr (pour l'arrêt des itérations
21
22*/
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
76    //if (Yt==TU+1)
77    // printf("%d\t%d\t%d\t%e\t%e (1)\n",Yk,Yj,Yi,YS1,YS2); 
78} //Yt>TU
79}
80
81//===========================================================================
82//                          methode backward
83
84backward (YREAL x1, YREAL x2, YREAL x3, YREAL x4, YREAL x5, YREAL x6, YREAL x7, YREAL x8) {
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}
94output 1 : gcx
95oupput 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      }*/
175}
Note: See TracBrowser for help on using the repository browser.