source: tags/version-1.1/trunk/src/solsor_dynspg_flt_.h @ 24

Last change on this file since 24 was 1, checked in by jbrlod, 13 years ago

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