source: trunk/src/solsor_gx2_x_r.h @ 49

/***************************************************************************
                 module classe solsor_gc1_x_r.h  -  description
***************************************************************************/
// Julien Brajard
// locean-ipsl.upmc, Paris, April 25, 2011
//===========================================================================
//                          methode forward
forward (YREAL x1, YREAL x2, YREAL x3, YREAL x4, YREAL x5, YREAL x6, YREAL x7, YREAL x8, YREAL x9)
{
  /*    1       from    gcx_dynspg_flt          1  i    j         t
        2       from    gcb_dynspg_flt          1  i    j         t
        3       from    solsor_gc1_x_r          1  i    j-1  k    t
        4       from    solsor_gc1_x_r          1  i-1  j    k    t
        5       from    solsor_gc1_x_r          1  i+1  j    k    t
        6       from    solsor_gc1_x_r          1  i    j+1  k    t
        7       from    solsor_gc1_x_r          1  i    j    k    t
        8       from    solsor_flag             2  NX   NY   k    t
        9       from    solsor_gc1_x_r          2  i    j    k-1  t
        } 
        output 1 : gcx
        oupput 2 : gcr (pour l'arrêt des itérations)
        
  */
  if(Yt==TU)
    {
      YS1=0;
      YS2=x9;
    }
  else {
    int t;
    if(Yt==TU+1 && neuler==0)
      t=0;
    else
      t=1;
    
    if (x8==1)
      {
        YS1=x7;
        YS2=x9;
      }
    else
      {
        if (Yk==0) {
          YS1=x7;
          YS2=x9;
        } //Yk=0;
        else {
          YREAL ztmp,zres;
          int ishift=(Yj-1)%2;
          if ((Yi+1)%2==ishift) { //! Guess black update
            ztmp = x2 
              - gcp(Yi,Yj,0,t) * x3   
              - gcp(Yi,Yj,1,t) * x4   
              - gcp(Yi,Yj,2,t) * x5   
              - gcp(Yi,Yj,3,t) * x6;
            //     ! Estimate of the residual
            zres = ztmp - x7;
            YS2 = zres * gcdmat(Yi,Yj,t) * zres;
            //               ! Guess update
            YS1 = sor * ztmp + (1.-sor) * x7;
          }
          else
            {
              YS2=x9;
              YS1=x7;
            }
          
          
          
        } // Yk>0
      } //if x8==0
    //if (Yt==TU+1) 
    //printf("%d\t%d\t%d\t%e\t%e (2)\n",Yk,Yj,Yi,YS1,YS2);  
  } //Yt>TU
}

//===========================================================================
//                          methode backward

backward (YREAL x1, YREAL x2, YREAL x3, YREAL x4, YREAL x5, YREAL x6, YREAL x7, YREAL x8, YREAL x9) {
  /*    1       from    gcx_dynspg_flt          1  i    j         t
        2       from    gcb_dynspg_flt          1  i    j         t
        3       from    solsor_gc1_x_r          1  i    j-1  k    t
        4       from    solsor_gc1_x_r          1  i-1  j    k    t
        5       from    solsor_gc1_x_r          1  i+1  j    k    t
        6       from    solsor_gc1_x_r          1  i    j+1  k    t
        7       from    solsor_gc1_x_r          1  i    j    k    t
        8       from    solsor_flag             2  NX   NY   k    t
        9       from    solsor_gc1_x_r          2  i    j    k-1  t
        } 
        output 1 : gcx
        oupput 2 : gcr (pour l'arrêt des itérations)
        
  */
  
  if(Yt==TU)
    {
      YJ2I9=0;
    }
  else {
    int t;
    if(Yt==TU+1 && neuler==0)
      t=0;
    else
      t=1;
    
    if (x8==1)
      {
        YJ1I7=1;
      }
    else
      {
        if (Yk==0) {
          YJ1I7=1;
        } //Yk=0;
        else {
          int ishift=(Yj-1)%2;
          if ((Yi+1)%2==ishift) { //! Guess black update
            YREAL dztmpdx2=1;
            YREAL dztmpdx3=-gcp(Yi,Yj,0,t);
            YREAL dztmpdx4=-gcp(Yi,Yj,1,t);
            YREAL dztmpdx5=-gcp(Yi,Yj,2,t);
            YREAL dztmpdx6=-gcp(Yi,Yj,3,t);
            /*ztmp = x2 
              - gcp(Yi,Yj,0,t) * x3   
              - gcp(Yi,Yj,1,t) * x4   
              - gcp(Yi,Yj,2,t) * x5   
              - gcp(Yi,Yj,3,t) * x6;*/
            
            YJ1I2=sor*dztmpdx2;
            YJ1I3=sor*dztmpdx3;
            YJ1I4=sor*dztmpdx4;
            YJ1I5=sor*dztmpdx5;
            YJ1I6=sor*dztmpdx6;
            YJ1I7=(1.-sor);          
            //       YS1 = sor * ztmp + (1.-sor) * x7;
          }
          else
            {
              YJ1I7=1;
              //  YS1=x7;
            }
        }          
      }
  }
  /*
  for (int i=0;i<  9;i++) 
    {
      
      for (int j=0;j<2;j++) {
        if (Yjac[j][i]>1 & YG1Y_solsor_gx2_x_r>1e-10) {
            printf("gx2 (%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]);
            printf(" YJ%dI%d=%e\n",j+1,i+1,Yjac[j][i]);
          }
      }
      }*/
}