/[escript]/branches/symbolic_from_3470/paso/src/FluxLimiter.c
ViewVC logotype

Diff of /branches/symbolic_from_3470/paso/src/FluxLimiter.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3867 by caltinay, Thu Feb 9 00:27:46 2012 UTC revision 3868 by caltinay, Thu Mar 15 06:07:08 2012 UTC
# Line 45  Paso_FCT_FluxLimiter* Paso_FCT_FluxLimit Line 45  Paso_FCT_FluxLimiter* Paso_FCT_FluxLimit
45          out->R=MEMALLOC(2*n,double);          out->R=MEMALLOC(2*n,double);
46          Esys_checkPtr(out->R);          Esys_checkPtr(out->R);
47    
         out->MQ_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),2*blockSize);  
48          out->R_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),2*blockSize);          out->R_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),2*blockSize);
49          out->u_tilde_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),blockSize);          out->u_tilde_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),blockSize);
50          out->antidiffusive_fluxes=Paso_SystemMatrix_alloc(fctp->transport_matrix->type,          out->antidiffusive_fluxes=Paso_SystemMatrix_alloc(fctp->transport_matrix->type,
# Line 53  Paso_FCT_FluxLimiter* Paso_FCT_FluxLimit Line 52  Paso_FCT_FluxLimiter* Paso_FCT_FluxLimit
52                              fctp->transport_matrix->row_block_size,                              fctp->transport_matrix->row_block_size,
53                              fctp->transport_matrix->col_block_size, TRUE);                              fctp->transport_matrix->col_block_size, TRUE);
54          out->borrowed_lumped_mass_matrix=fctp->lumped_mass_matrix;          out->borrowed_lumped_mass_matrix=fctp->lumped_mass_matrix;
           
55      }      }
56      if (Esys_noError()) {      if (Esys_noError()) {
57          return out;          return out;
# Line 71  void Paso_FCT_FluxLimiter_free(Paso_FCT_ Line 69  void Paso_FCT_FluxLimiter_free(Paso_FCT_
69         MEMFREE(in->u_tilde);         MEMFREE(in->u_tilde);
70         MEMFREE(in->MQ);         MEMFREE(in->MQ);
71         MEMFREE(in->R);         MEMFREE(in->R);
        Paso_Coupler_free(in->MQ_coupler);  
72         Paso_Coupler_free(in->R_coupler);         Paso_Coupler_free(in->R_coupler);
73         Paso_Coupler_free(in->u_tilde_coupler);         Paso_Coupler_free(in->u_tilde_coupler);
74         MEMFREE(in);         MEMFREE(in);
# Line 88  void Paso_FCT_FluxLimiter_setU_tilda(Pas Line 85  void Paso_FCT_FluxLimiter_setU_tilda(Pas
85    const double *lumped_mass_matrix = flux_limiter->borrowed_lumped_mass_matrix;    const double *lumped_mass_matrix = flux_limiter->borrowed_lumped_mass_matrix;
86    index_t i, iptr_ij;    index_t i, iptr_ij;
87    
88    #pragma omp parallel private(i, iptr_ij)    #pragma omp parallel private(i)
89    {    {
90      #pragma omp for schedule(static)      #pragma omp for schedule(static)
91      for (i = 0; i < n; ++i) {      for (i = 0; i < n; ++i) {
92             const double m=lumped_mass_matrix[i];                  const double m=lumped_mass_matrix[i];
93             flux_limiter->u_tilde[i]=Mu_tilda[i]/m;          if (m > 0 ) {
94                        flux_limiter->u_tilde[i]=Mu_tilda[i]/m;
95            } else {
96                flux_limiter->u_tilde[i]=Mu_tilda[i];
97            }
98          
99      }              }        
100          }
101      /* distribute u_tilde: */    /* distribute u_tilde: */
102      #pragma omp master    Paso_Coupler_startCollect(flux_limiter->u_tilde_coupler,flux_limiter->u_tilde);
103      {    /*
         Paso_Coupler_startCollect(flux_limiter->u_tilde_coupler,flux_limiter->u_tilde);  
     }  
     /*  
104       * calculate MQ_P[i] = lumped_mass_matrix[i] * max_{j} (\tilde{u}[j]- \tilde{u}[i] )       * calculate MQ_P[i] = lumped_mass_matrix[i] * max_{j} (\tilde{u}[j]- \tilde{u}[i] )
105       *           MQ_N[i] = lumped_mass_matrix[i] * min_{j} (\tilde{u}[j]- \tilde{u}[i] )       *           MQ_N[i] = lumped_mass_matrix[i] * min_{j} (\tilde{u}[j]- \tilde{u}[i] )
106       *       *
107       */       */
108      /* first we calculate the min and max of u_tilda in the main block    /* first we calculate the min and max of u_tilda in the main block
109         QP, QN are used to hold the result */         QP, QN are used to hold the result */
110          #pragma omp parallel private(i)
111       #pragma omp  for schedule(static)    {
      for (i = 0; i < n; ++i) {  
         double u_min_i=LARGE_POSITIVE_FLOAT;  
         double u_max_i=-LARGE_POSITIVE_FLOAT;  
         #pragma ivdep  
         for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {  
              const index_t j=pattern->mainPattern->index[iptr_ij];  
              const double u_j=flux_limiter->u_tilde[j];  
              u_min_i=MIN(u_min_i,u_j);  
              u_max_i=MAX(u_max_i,u_j);  
         }  
         flux_limiter->MQ[2*i] = u_min_i;    
         flux_limiter->MQ[2*i+1] = u_max_i;  
      }    
      /* complete distribute u_tilde: */  
      #pragma omp master  
      {    
         Paso_Coupler_finishCollect(flux_limiter->u_tilde_coupler);  
     remote_u_tilde=Paso_Coupler_borrowRemoteData(flux_limiter->u_tilde_coupler);  
      }  
      /* now we look at the couple matrix and set the final value for QP, QN */  
112       #pragma omp  for schedule(static)       #pragma omp  for schedule(static)
113       for (i = 0; i < n; ++i) {       for (i = 0; i < n; ++i) {
114          const double u_i = flux_limiter->u_tilde[i];          if (  flux_limiter->borrowed_lumped_mass_matrix[i]> 0) { /* no constraint */    
115          double u_min_i = flux_limiter->MQ[2*i];          double u_min_i=LARGE_POSITIVE_FLOAT;
116          double u_max_i = flux_limiter->MQ[2*i+1];          double u_max_i=-LARGE_POSITIVE_FLOAT;
117          #pragma ivdep          #pragma ivdep
118          for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
119               const index_t j  = pattern->col_couplePattern->index[iptr_ij];              const index_t j=pattern->mainPattern->index[iptr_ij];
120               const double u_j = remote_u_tilde[j];              const double u_j=flux_limiter->u_tilde[j];
121               u_min_i=MIN(u_min_i,u_j);              u_min_i=MIN(u_min_i,u_j);
122               u_max_i=MAX(u_max_i,u_j);              u_max_i=MAX(u_max_i,u_j);
123           }          }
124           flux_limiter->MQ[2*i  ] = ( u_min_i-u_i ) * lumped_mass_matrix[i];  /* M_C*Q_min*/          flux_limiter->MQ[2*i] = u_min_i;  
125           flux_limiter->MQ[2*i+1] = ( u_max_i-u_i ) * lumped_mass_matrix[i] ;  /* M_C*Q_max */          flux_limiter->MQ[2*i+1] = u_max_i;
126            
127        } else {
128                flux_limiter->MQ[2*i  ] =LARGE_POSITIVE_FLOAT;
129                flux_limiter->MQ[2*i+1] =LARGE_POSITIVE_FLOAT;
130        }
131       }       }
132    } /* parallel region */    } /* parallel region */
133      /* complete distribute u_tilde: */
134      Paso_Coupler_finishCollect(flux_limiter->u_tilde_coupler);
135      remote_u_tilde=Paso_Coupler_borrowRemoteData(flux_limiter->u_tilde_coupler);
136    
137      /* now we look at the couple matrix and set the final value for QP, QN */
138      #pragma omp parallel private(i, iptr_ij)
139      {
140          #pragma omp  for schedule(static)
141          for (i = 0; i < n; ++i) {
142            if ( flux_limiter->borrowed_lumped_mass_matrix[i]> 0) { /* no constraint */
143            const double u_i = flux_limiter->u_tilde[i];
144            double u_min_i = flux_limiter->MQ[2*i];
145            double u_max_i = flux_limiter->MQ[2*i+1];
146            #pragma ivdep
147            for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
148                const index_t j  = pattern->col_couplePattern->index[iptr_ij];
149                const double u_j = remote_u_tilde[j];
150                u_min_i=MIN(u_min_i,u_j);
151                u_max_i=MAX(u_max_i,u_j);
152            }
153            flux_limiter->MQ[2*i  ] = ( u_min_i-u_i ) * lumped_mass_matrix[i];  /* M_C*Q_min*/
154            flux_limiter->MQ[2*i+1] = ( u_max_i-u_i ) * lumped_mass_matrix[i] ;  /* M_C*Q_max */
155            }
156    
157           }
158    Paso_Coupler_startCollect(flux_limiter->MQ_coupler,flux_limiter->MQ);    } /* parallel region */
   Paso_Coupler_finishCollect(flux_limiter->MQ_coupler);  
159  }  }
160    
161  /* starts to update a vector (not given yet) from the antidiffusion fluxes in flux_limiter->antidiffusive_fluxes  /* starts to update a vector (not given yet) from the antidiffusion fluxes in flux_limiter->antidiffusive_fluxes
# Line 162  void Paso_FCT_FluxLimiter_addLimitedFlux Line 167  void Paso_FCT_FluxLimiter_addLimitedFlux
167    const double* u_tilde = flux_limiter->u_tilde;    const double* u_tilde = flux_limiter->u_tilde;
168    const double* remote_u_tilde=Paso_Coupler_borrowRemoteData(flux_limiter->u_tilde_coupler);    const double* remote_u_tilde=Paso_Coupler_borrowRemoteData(flux_limiter->u_tilde_coupler);
169    Paso_SystemMatrix * adf=flux_limiter->antidiffusive_fluxes;    Paso_SystemMatrix * adf=flux_limiter->antidiffusive_fluxes;
   /* const double *lumped_mass_matrix = flux_limiter->borrowed_lumped_mass_matrix; */  
170    index_t iptr_ij;    index_t iptr_ij;
171    dim_t i;    dim_t i;
172    
# Line 170  void Paso_FCT_FluxLimiter_addLimitedFlux Line 174  void Paso_FCT_FluxLimiter_addLimitedFlux
174    {    {
175         #pragma omp for schedule(static)           #pragma omp for schedule(static)  
176         for (i = 0; i < n; ++i) {         for (i = 0; i < n; ++i) {
       
           const double u_tilde_i=u_tilde[i];  
           double P_P_i=0.;  
           double P_N_i=0.;  
       const double MQ_min=flux_limiter->MQ[2*i];  
       const double MQ_max=flux_limiter->MQ[2*i+1];  
177        double R_N_i =1;        double R_N_i =1;
178        double R_P_i =1;        double R_P_i =1;
179            #pragma ivdep        if ( flux_limiter->borrowed_lumped_mass_matrix[i]> 0) { /* no constraint */
180            for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {          const double u_tilde_i=u_tilde[i];
181               const index_t j=pattern->mainPattern->index[iptr_ij];          double P_P_i=0.;
182           if (i != j ) {          double P_N_i=0.;
183                      const double f_ij=adf->mainBlock->val[iptr_ij];          const double MQ_min=flux_limiter->MQ[2*i];
184                  const double u_tilde_j=u_tilde[j];          const double MQ_max=flux_limiter->MQ[2*i+1];
185                  /* pre-limiter */          #pragma ivdep
186                      if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
187                         adf->mainBlock->val[iptr_ij]=0;            const index_t j=pattern->mainPattern->index[iptr_ij];
188                  } else {            if (i != j ) {
189                          if (f_ij <=0) {                const double f_ij=adf->mainBlock->val[iptr_ij];
190                              P_N_i+=f_ij;                const double u_tilde_j=u_tilde[j];
191                          } else {                /* pre-limiter */
192                              P_P_i+=f_ij;                if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {
193                          }                  adf->mainBlock->val[iptr_ij]=0;
194                      }                } else {
195           }                    if (f_ij <=0) {
196        }                    P_N_i+=f_ij;
197                            } else {
198        /* now the couple matrix: */                    P_P_i+=f_ij;
199        #pragma ivdep                    }
200        for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {                }
201               const index_t j=pattern->col_couplePattern->index[iptr_ij];            }
202               const double f_ij=adf->col_coupleBlock->val[iptr_ij];          }
203           const double u_tilde_j=remote_u_tilde[j];          
204            /* pre-limiter */          /* now the couple matrix: */
205                if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {          #pragma ivdep
206                     adf->col_coupleBlock->val[iptr_ij]=0;          for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
207            } else {            const index_t j=pattern->col_couplePattern->index[iptr_ij];
208                    if (f_ij <=0) {            const double f_ij=adf->col_coupleBlock->val[iptr_ij];
209                        P_N_i+=f_ij;            const double u_tilde_j=remote_u_tilde[j];
210                    } else {              /* pre-limiter */
211                        P_P_i+=f_ij;              if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {
212                    }              adf->col_coupleBlock->val[iptr_ij]=0;
213                }              } else {
214        }              if (f_ij <=0) {
215        /* finally the R+ and R- are calculated */                  P_N_i+=f_ij;
216                } else {
217                    P_P_i+=f_ij;
218                }
219                }
220            }
221            /* finally the R+ and R- are calculated */
222    
223            if (P_N_i<0) R_N_i=MIN(1,MQ_min/P_N_i);          if (P_N_i<0) R_N_i=MIN(1,MQ_min/P_N_i);
224            if (P_P_i>0) R_P_i=MIN(1,MQ_max/P_P_i);          if (P_P_i>0) R_P_i=MIN(1,MQ_max/P_P_i);
225          }
226            flux_limiter->R[2*i]   = R_N_i;            flux_limiter->R[2*i]   = R_N_i;
227        flux_limiter->R[2*i+1] = R_P_i;        flux_limiter->R[2*i+1] = R_P_i;
228        }        }

Legend:
Removed from v.3867  
changed lines
  Added in v.3868

  ViewVC Help
Powered by ViewVC 1.1.26