source: trunk/src/solsor_dynspg_flt.h @ 54

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

option SOLOSORYAO off
  • Property svn:eol-style set to native
File size: 9.2 KB
Line 
1
2/***************************************************************************
3                 module classe solsor_dynspg_flt.h  -  description
4***************************************************************************/
5// Mohamed Berrada
6// locean-ipsl.upmc, Paris, April 23, 2009
7//===========================================================================
8//                          methode forward
9forward ()
10{
11  // FILE * fid;
12 
13  if(Yt==TU)
14    return;
15  else{
16    if(Yi==0 && Yj==0){
17      int t,jn,jj,ji, icount = 0,ncut=0,niter=0,ishift;
18      niter=0;ncut=0;
19      if(Yt==TU+1 && neuler==0)
20        t=0;
21      else
22        t=1;
23      //---------------------
24      //initialisation de gcx
25     
26      /*test solsor
27      if (Yt==TU+1) {
28       
29        fid=fopen("../data_out/test_solsor_mb","w");
30        }*/
31     
32     
33      for(jj=0;jj<NY;jj++)
34        for(ji=0;ji<NX;ji++){
35          //    gcxb(ji,jj)=gcx(ji,jj);
36          gcx(ji,jj)=YS_gcx_dynspg_flt(0,ji,jj,Yt);
37 
38          /*test solsor
39          if (Yt==TU+1) {
40           
41            fprintf(fid,"%d\t%d\t%d\t%e\t%e\n",0,jj,ji,gcx(ji,jj),gcr(ji,jj));
42            }*/
43         
44
45        }
46      //if(Yt==6) printf("------------------gcx(ji,jj,0,t)  =%lf\n",gcx(28,2) );
47      //---------------------
48      //boundarie conditions -- ceci doit etre a l'interieur de loop nmax mais cas special 
49      for(jj=0;jj<NY;jj++){// ! East-West boundaries
50        gcx(0,jj)=0.;gcx(NX-1,jj)=0.;
51      }
52      for(ji=0;ji<NX;ji++){// ! North-South boundaries
53        gcx(ji,0)=0.;gcx(ji,NY-1)=0.;
54      }
55     
56     
57      double ztmp,zres,zres2,res=0.;
58      //---------------------
59      //boucle
60      //xsolmat_init();
61      for(jn=1;jn<=nmax;jn++){
62        // cas special        ! applied the lateral boundary conditions
63        // CALL lbc_lnk_e( gcx, c_solver_pt, 1. ) ;  //c_solver_pt='T'   
64        //         ! Residus   
65        //         ! Guess black update
66       
67        for(jj=1;jj<NY-1;jj++){
68          ishift = jj%2;
69          for(ji=1+ishift;ji<NX-1;ji+=2){
70            ztmp = YS_gcb_dynspg_flt(0,ji,jj,Yt) 
71              - gcp(ji,jj,0,t) * gcx(ji  ,jj-1)   
72              - gcp(ji,jj,1,t) * gcx(ji-1,jj  )   
73              - gcp(ji,jj,2,t) * gcx(ji+1,jj  )   
74              - gcp(ji,jj,3,t) * gcx(ji  ,jj+1);
75            //     ! Estimate of the residual
76            zres = ztmp - gcx(ji,jj);
77            gcr(ji,jj) = zres * gcdmat(ji,jj,t) * zres;
78            //               ! Guess update
79            gcx(ji,jj) = sor * ztmp + (1.-sor) * gcx(ji,jj);
80 
81
82            /*test solsor
83            if (Yt==TU+1) {
84             
85              fprintf(fid,"%d\t%d\t%d\t%e\t%e\n",jn,jj,ji,gcx(ji,jj),gcr(ji,jj));
86            }
87            */
88
89          }
90        }
91        icount = icount + 1; 
92       
93        //cas special     ! applied the lateral boundary conditions
94        //   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )  ;
95       
96        //   ! Guess red update
97        for(jj=1;jj<NY-1;jj++){
98          ishift = (jj-1)%2;
99          for(ji=1+ishift;ji<NX-1;ji+=2){
100            ztmp = YS_gcb_dynspg_flt(0,ji,jj,Yt) 
101              - gcp(ji,jj,0,t) * gcx(ji  ,jj-1)   
102              - gcp(ji,jj,1,t) * gcx(ji-1,jj  )   
103              - gcp(ji,jj,2,t) * gcx(ji+1,jj  )   
104              - gcp(ji,jj,3,t) * gcx(ji  ,jj+1);
105            //     ! Estimate of the residual
106            zres = ztmp - gcx(ji,jj);
107            gcr(ji,jj) = zres * gcdmat(ji,jj,t) * zres;
108            //               ! Guess update
109            gcx(ji,jj) = sor * ztmp + (1.-sor) * gcx(ji,jj);
110
111            /*test solsor
112            if (Yt==TU+1) {
113             
114              fprintf(fid,"%d\t%d\t%d\t%e\t%e\n",jn,jj,ji,gcx(ji,jj),gcr(ji,jj));
115              }*/
116          }
117        }
118        //      if(jn>=nmax-1) xtest("gcx__",gcx(28,3),28,3,0,Yt);
119        icount = icount + 1;
120        //     ! test of convergence
121        if( jn > nmin && (jn-nmin)%nmod == 0 ) {
122          zres2=0.;
123          for(jj=1;jj<NY-1;jj++)
124            for(ji=1;ji<NX-1;ji++)
125              if(zres2<gcr(ji,jj)) zres2=gcr(ji,jj);
126          //      zres2 = MAXVAL( gcr(2:nlci-1,2:nlcj-1) );
127          //               ! test of convergence
128          if( zres2 < resmax || jn == nmax ){
129            res = sqrt( zres2 );
130            niter = jn; tniter[Yt]=jn;
131            ncut = 999;
132            //        printf("FOR=======tniter(%d)=%d\n",Yt,tniter[Yt]);
133          }
134        }
135        //      if(jn>=nmax-1) xtest("gcx____",gcx(28,3),28,3,0,Yt);
136       
137        //         ! indicator of non-convergence or explosion
138        //         IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
139        if( ncut == 999 ) break;
140      }
141      //    cout<<"niter="<<niter<<endl;
142      //    printf("res=%24.16e\n",res);
143      // CALL lbc_lnk_e( gcx, c_solver_pt, 1. );//    ! boundary conditions
144      /* FILE* fp;
145         fp=fopen("testh.pof","w");
146         for(int i=0;i<NX;i++)
147         for(int j=0;j<NY;j++){
148         fprintf(fp,"%4d  %4d %4d  %24.16e\n",Yt+1,i,j,gcr(i,j));
149         }
150         fclose(fp);*/
151      /*      xtest("gcx",gcx(28,3),28,3,0,Yt);
152              xtest("gcb",YS_gcb_dynspg_flt(0,28,3,Yt),28,3,0,Yt); */
153    /*test solsor
154    if (Yt==TU+1) {
155      fclose(fid);
156    }
157    */
158   
159    }
160  }
161  //
162}
163
164//===========================================================================
165//                          methode backward
166
167backward ()
168{
169
170  if(Ycurward==BACKWARD){
171    if(Yt==TU)
172      return;
173    else{
174      if(Yi==0 && Yj==0){
175        double G_ztmp;
176        int t,jn,jj,ji, ishift,istart;
177        if(Yt==TU+1 && neuler==0)
178          t=0;
179        else
180          t=1;
181        //---------------------
182        for(ji=0;ji<NX;ji++){
183          for(jj=0;jj<NY;jj++){
184            YG_gcx_dynspg_flt(0,ji,jj,Yt)=0.;
185            YG_gcb_dynspg_flt(0,ji,jj,Yt)=0.;
186          }
187        }
188        for(jj=0;jj<NY;jj++){
189          G_gcx(0,jj)=0.;G_gcx(NX-1,jj)=0.;
190          G_gcx0(0,jj)=0.;G_gcx0(NX-1,jj)=0.;
191        }
192        for(ji=0;ji<NX;ji++){
193          G_gcx(ji,0)=0.;G_gcx(ji,NY-1)=0.;
194          G_gcx0(ji,0)=0.;G_gcx0(ji,NY-1)=0.;
195        }
196
197        for(ji=1;ji<NX-1;ji++)
198          for(jj=1;jj<NY-1;jj++){
199            //      YG_gcx2(0,ji,jj,Yt)=Ytb_gcx(ji,jj);
200            G_gcx(ji,jj)=YG_gcx2(0,ji,jj,Yt);
201            G_gcx0(ji,jj)=0.;
202          }
203
204        for(jn=tniter[Yt];jn>=1;jn--){
205          //    printf("BACK1=======YG_gcx_dynspg_flt(%d)=%23.16e\n",Yt,G_gcx(9,9));
206          //    printf("BACK2=======YG_gcb_dynspg_flt(%d)=%23.16e\n",Yt,YG_gcb_dynspg_flt(0,9,9,Yt));
207
208          //  BACK   ! Guess red update
209          for(jj=NY-2;jj>=1;jj--){
210            ishift = (jj-1)%2;
211            istart=NX-2-(ishift-1)%2*(NX%2-1)%2-ishift*NX%2;
212            for(ji=istart;ji>=1;ji-=2){
213              G_ztmp=G_gcx(ji,jj)*sor;
214              G_gcx0(ji,jj)=G_gcx(ji,jj)*(1.-sor);
215              YG_gcb_dynspg_flt(0,ji,jj,Yt)+=G_ztmp;
216              G_gcx0(ji  ,jj-1)+=- G_ztmp*gcp(ji,jj,0,t) ;
217              G_gcx0(ji-1,jj  )+=- G_ztmp*gcp(ji,jj,1,t) ;
218              G_gcx0(ji+1,jj  )+=- G_ztmp*gcp(ji,jj,2,t) ;
219              G_gcx0(ji  ,jj+1)+=- G_ztmp*gcp(ji,jj,3,t) ;
220            }
221          }
222
223          //   BACK      ! Guess black update
224         
225          for(jj=NY-2;jj>=1;jj--){
226            ishift = jj%2;
227            istart=NX-2-(ishift-1)%2*(NX%2-1)%2-ishift*NX%2;
228            for(ji=istart;ji>=1;ji-=2){
229              G_ztmp=(G_gcx(ji,jj)+G_gcx0(ji,jj))*sor;
230              G_gcx0(ji,jj)=(G_gcx(ji,jj)+G_gcx0(ji,jj))*(1.-sor);
231              YG_gcb_dynspg_flt(0,ji,jj,Yt)+=G_ztmp;
232              G_gcx0(ji  ,jj-1)+=- G_ztmp*gcp(ji,jj,0,t) ;
233              G_gcx0(ji-1,jj  )+=- G_ztmp*gcp(ji,jj,1,t) ;
234              G_gcx0(ji+1,jj  )+=- G_ztmp*gcp(ji,jj,2,t) ;
235              G_gcx0(ji  ,jj+1)+=-G_ztmp* gcp(ji,jj,3,t) ;
236            }
237          }
238          for(jj=0;jj<NY;jj++)
239            for(ji=0;ji<NX;ji++){
240              G_gcx(ji,jj)=G_gcx0(ji,jj);
241              G_gcx0(ji,jj)=0.;
242            }
243          for(jj=0;jj<NY;jj++){
244            G_gcx0(0,jj)=0.;G_gcx0(NX-1,jj)=0.;
245          }
246          for(ji=0;ji<NX;ji++){
247            G_gcx0(ji,0)=0.;G_gcx0(ji,NY-1)=0.;
248          }
249         
250         
251        }
252        for(jj=1;jj<NY-1;jj++)
253          for(ji=1;ji<NX-1;ji++){
254            YG_gcx_dynspg_flt(0,ji,jj,Yt)+=G_gcx(ji,jj);
255          }
256       
257
258       
259      }
260    }
261   
262  }
263  else if(Ycurward==LINWARD){
264    if(Yt==TU)
265      return;
266    else{
267      if(Yi==0){
268        int t,jn,jj,ji, ishift;
269        if(Yt==TU+1 && neuler==0)
270          t=0;
271        else
272          t=1;
273        //---------------------
274        //initialisation de gcx
275        for(jj=0;jj<NY;jj++)
276          for(ji=0;ji<NX;ji++){
277            Ytb_gcx(ji,jj)=YG_gcx_dynspg_flt(0,ji,jj,Yt);
278          }
279        //---------------------
280        //boundarie conditions -- ceci doit etre a l'interieur de loop nmax mais cas special 
281        for(jj=0;jj<NY;jj++){// ! East-West boundaries
282          Ytb_gcx(0,jj)=0.;Ytb_gcx(NX-1,jj)=0.;
283        }
284        for(ji=0;ji<NX;ji++){// ! North-South boundaries
285          Ytb_gcx(ji,0)=0.;Ytb_gcx(ji,NY-1)=0.;
286        }
287       
288        double Ytb_ztmp;
289        //---------------------
290        //boucle
291        //xsolmat_init();
292        //      printf("LIN=======tniter(%d)=%d\n",Yt,tniter[Yt]);
293        for(jn=1;jn<=tniter[Yt];jn++){
294          // cas special        ! applied the lateral boundary conditions
295          // CALL lbc_lnk_e( gcx, c_solver_pt, 1. ) ;  //c_solver_pt='T'   
296          //         ! Residus   
297          //         ! Guess black update
298         
299          for(jj=1;jj<NY-1;jj++){
300            ishift = jj%2;//MOD( jj, 2 );
301            for(ji=1+ishift;ji<NX-1;ji+=2){
302              Ytb_ztmp = YG_gcb_dynspg_flt(0,ji,jj,Yt) 
303                - gcp(ji,jj,0,t) * Ytb_gcx(ji  ,jj-1)   
304                - gcp(ji,jj,1,t) * Ytb_gcx(ji-1,jj  )   
305                - gcp(ji,jj,2,t) * Ytb_gcx(ji+1,jj  )   
306                - gcp(ji,jj,3,t) * Ytb_gcx(ji  ,jj+1);
307              //               ! Guess update
308              Ytb_gcx(ji,jj) = sor * Ytb_ztmp + (1.-sor) * Ytb_gcx(ji,jj);
309            }
310          }
311         
312          //cas special     ! applied the lateral boundary conditions
313          //   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )  ;
314         
315          //   ! Guess red update
316          for(jj=1;jj<NY-1;jj++){
317            ishift = (jj-1)%2;
318            for(ji=1+ishift;ji<NX-1;ji+=2){
319              Ytb_ztmp = YG_gcb_dynspg_flt(0,ji,jj,Yt) 
320                - gcp(ji,jj,0,t) * Ytb_gcx(ji  ,jj-1)   
321                - gcp(ji,jj,1,t) * Ytb_gcx(ji-1,jj  )   
322                - gcp(ji,jj,2,t) * Ytb_gcx(ji+1,jj  )   
323                - gcp(ji,jj,3,t) * Ytb_gcx(ji  ,jj+1);
324              //               ! Guess update
325              Ytb_gcx(ji,jj) = sor * Ytb_ztmp + (1.-sor) * Ytb_gcx(ji,jj);
326            }
327          }
328         
329        }
330        for(ji=1;ji<NX-1;ji++)
331          for(jj=1;jj<NY-1;jj++)
332            YG_gcx2(0,ji,jj,Yt)=Ytb_gcx(ji,jj);
333      }
334    }
335    //
336    return;
337  }
338  //
339}
340
341//===========================================================================
342//**********************  FIN DU MODULE solsor_dynspg_flt   *****************
Note: See TracBrowser for help on using the repository browser.