/[escript]/trunk/paso/src/SolverFCT_FluxControl.c
ViewVC logotype

Diff of /trunk/paso/src/SolverFCT_FluxControl.c

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

trunk/paso/src/Solver_FluxControl.c revision 1361 by gross, Fri Dec 14 09:26:51 2007 UTC trunk/paso/src/SolverFCT_FluxControl.c revision 1408 by gross, Mon Feb 4 07:19:50 2008 UTC
# Line 26  Line 26 
26  #include "SolverFCT.h"  #include "SolverFCT.h"
27  #include "PasoUtil.h"  #include "PasoUtil.h"
28    
 #define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.)  
 #define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b)))  
 #define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b))))  
   
 #define FLUX_L(a,b) SUPERBEE(a,b)  /* alter for other flux limiter */  
   
 #define FLUX_LIMITER(a) FLUX_L(a,1)  
29    
30  /**************************************************************/  /**************************************************************/
31    
32  /* free all memory used by FluxControl                                */  /* create the low order transport matrix and stores it negative values
33     * into the iteration_matrix accept the main diagonal which is stored seperately
34  void Paso_Solver_FluxControl_free(Paso_Solver_FluxControl* in) {   * if fc->iteration_matrix==NULL, fc->iteration_matrix is allocated
35       if (in!=NULL) {   *
36          Paso_SystemMatrix_freeBuffer(in->matrix);   * a=transport_matrix
37          Paso_SystemMatrix_free(in->matrix);   * b= low_order_transport_matrix = - iteration_matrix
38          MEMFREE(in->colorOf);   * c=main diagonal low_order_transport_matrix
39          MEMFREE(in->main_iptr);   * initialize c[i] mit a[i,i]
40          MEMFREE(in);   *    d_ij=max(0,-a[i,j],-a[j,i])  
41       }   *    b[i,j]=-(a[i,j]+d_ij)        
42  }   *    c[i]-=d_ij                  
43     */
44    
45  /**************************************************************/  void Paso_FCTransportProblem_setLowOrderOperator(Paso_FCTransportProblem * fc) {
46      dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i;
47  /*   constructs a flux control mechanism                      */    index_t color, iptr_ij,j,iptr_ji;
48      Paso_SystemMatrixPattern *pattern;
49  Paso_Solver_FluxControl* Paso_SolverFCT_getFluxControl(Paso_SystemMatrix * A) {    register double d_ij, sum, rtmp1, rtmp2;
50    
51    Paso_Solver_FluxControl* out=NULL;    if (fc==NULL) return;
52    dim_t n,i;    if (fc->iteration_matrix==NULL) {
53    index_t iptr,iptr_main,k;        fc->iteration_matrix=Paso_SystemMatrix_alloc(fc->transport_matrix->type,
54                                                       fc->transport_matrix->pattern,
55    if (A==NULL) return out;                                                     fc->transport_matrix->row_block_size,
56    n=Paso_SystemMatrix_getTotalNumRows(A);                                                     fc->transport_matrix->col_block_size);
57    if (A->block_size!=1) {    }
         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_getFluxControl: block size > 1 is not supported.");  
         return NULL;  
   }  
   out=MEMALLOC(1,Paso_Solver_FluxControl);  
   if (Paso_checkPtr(out)) return NULL;  
   
   out->matrix=Paso_SystemMatrix_reference(A);  
   out->colorOf=NULL;  
   out->main_iptr=NULL;  
     
   
   /* allocations: */    
   out->colorOf=MEMALLOC(n,index_t);  
   out->main_iptr=MEMALLOC(n,index_t);  
   if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) ) ) {  
       printf("Paso_SolverFCT_getFluxControl: Revise coloring!!\n");  
       Paso_Pattern_color(A->mainBlock->pattern,&(out->num_colors),out->colorOf);  
       Paso_SystemMatrix_allocBuffer(A);  
58    
59        #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)    if (Paso_noError()) {
60          pattern=fc->iteration_matrix->pattern;
61          n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
62          #pragma omp parallel for private(i,sum,iptr_ij,j,iptr_ji,rtmp1, rtmp2,d_ij)  schedule(static)
63        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
64          for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
65               iptr_main=A->mainBlock->pattern->ptr[0]-1;            /* look at a[i,j] */
66                for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; iptr++) {            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
67                     if (A->mainBlock->pattern->index[iptr]==i) {                j=pattern->mainPattern->index[iptr_ij];
68                          iptr_main=iptr;                if (j!=i) {
69                     /* find entry a[j,i] */
70                     for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
71                        if (pattern->mainPattern->index[iptr_ji]==i) {
72                            rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];
73                            rtmp2=fc->transport_matrix->mainBlock->val[iptr_ji];
74                            d_ij=-MIN3(0.,rtmp1,rtmp2);
75                            fc->iteration_matrix->mainBlock->val[iptr_ij]=-(rtmp1+d_ij);
76                            sum-=d_ij;
77                          break;                          break;
78                     }                      }
79                 }                   }
80                 out->main_iptr[i]=iptr_main;               }
81                 if (iptr_main==A->mainBlock->pattern->ptr[0]-1)           }
82                    Paso_setError(VALUE_ERROR, "Paso_SolverFCT_getFluxControl: no main diagonal");           /* TODO process couple block */
83             }          
84         }           /* set main diagonal entry */
85             fc->main_diagonal_low_order_transport_matrix[i]=sum;
86    }        }
   if (Paso_noError()) {  
      return out;  
   } else {  
      Paso_Solver_FluxControl_free(out);  
      return NULL;  
87    }    }
88  }  }
89    /*
90     * out_i=m_i u_i + alpha \sum_{j <> i} l_{ij} (u_j-u_i) + beta q_i
91     *
92     */
93    void Paso_SolverFCT_setMuPaLuPbQ(double* out,
94                                     const double* M,
95                                     const double* u,
96                                     const double a,
97                                     Paso_SystemMatrix *L,
98                                     const double b,
99                                     const double* Q)
100    {
101      dim_t n, i, j;
102      Paso_SystemMatrixPattern *pattern;
103      double *remote_u;
104      register double sum, u_i, l_ij;
105      register index_t iptr_ij;
106    
107  /**************************************************************/    n=Paso_SystemMatrix_getTotalNumRows(L);
108    
109  /* adds A plus stabelising diffusion into the matrix B        */    if (ABS(a)>0) {
110          Paso_SystemMatrix_startCollect(L,u);
111      }
112      #pragma omp parallel for private(i) schedule(static)
113      for (i = 0; i < n; ++i) {
114          out[i]=M[i]*u[i]+b*Q[i];
115      }
116      if (ABS(a)>0) {
117          pattern=L->pattern;
118          remote_u=Paso_SystemMatrix_finishCollect(L);
119          #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
120          for (i = 0; i < n; ++i) {
121              sum=0;
122              u_i=u[i];
123              #pragma ivdep
124          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
125                   j=pattern->mainPattern->index[iptr_ij];
126                   l_ij=L->mainBlock->val[iptr_ij];
127                   sum+=l_ij*(u[j]-u_i);
128              }
129              #pragma ivdep
130          for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
131                   j=pattern->couplePattern->index[iptr_ij];
132                   l_ij=L->coupleBlock->val[iptr_ij];
133                   sum+=l_ij*(remote_u[j]-u_i);
134              }
135              out[i]+=a*sum;
136          }
137      }
138    }
139    /*
140     * calculate QP[i] max_{i in L->pattern[i]} (u[j]-u[i] )
141     *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
142     *
143    */
144    void Paso_SolverFCT_setQs(const double* u,double* QN, double* QP, Paso_SystemMatrix *L)
145    {
146      dim_t n, i, j;
147      Paso_SystemMatrixPattern *pattern;
148      double *remote_u;
149      register double u_i, u_min_i, u_max_i, u_j;
150      register index_t iptr_ij;
151    
152      Paso_SystemMatrix_startCollect(L,u);
153      n=Paso_SystemMatrix_getTotalNumRows(L);
154      pattern=L->pattern;
155      remote_u=Paso_SystemMatrix_finishCollect(L);
156      #pragma omp parallel for schedule(static) private(i, u_i, u_min_i, u_max_i, j, u_j, iptr_ij)
157      for (i = 0; i < n; ++i) {
158         u_i=u[i];
159         u_min_i=u_i;
160         u_max_i=u_i;
161         #pragma ivdep
162         for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
163             j=pattern->mainPattern->index[iptr_ij];
164             u_j=u[j];
165             u_min_i=MIN(u_min_i,u_j);
166             u_max_i=MIN(u_max_i,u_j);
167         }
168         #pragma ivdep
169         for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
170              j=pattern->couplePattern->index[iptr_ij];
171              u_j=remote_u[j];
172              u_min_i=MIN(u_min_i,u_j);
173              u_max_i=MIN(u_max_i,u_j);
174          }
175          QN[i]=u_min_i-u_i;
176          QP[i]=u_max_i-u_i;
177      }
178    }
179    
180  /* d_ij=alpha*max(0,-a[i,j],-a[j,i])  */  /*
181  /* b[i,j]+=alpha*(a[i,j]+d_ij)  */   *
182  /* b[j,i]+=alpha*(a[j,i]+d_ij)  */   * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])
183  /* b[i,i]-=alpha*d_ij  */   *
184  /* b[j,j]-=alpha*d_ij  */   * m=fc->mass matrix
185     * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
186     */
187    void Paso_FCTransportProblem_updateAntiDiffusionFlux(const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix,const double a, const double b, const double* u)
188    {
189      dim_t n, j, i;
190      double *remote_u;
191      index_t iptr_ij;
192      const double b2=-b;
193      register double m_ij, ml_ij, k_ij, u_i;
194      Paso_SystemMatrixPattern *pattern;
195      Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);
196      n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
197      pattern=fc->iteration_matrix->pattern;
198      remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
199    
200      if ( (ABS(a) >0 ) ) {
201          if ( (ABS(b) >0 ) ) {
202             #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij)
203             for (i = 0; i < n; ++i) {
204                u_i=u[i];
205                #pragma ivdep
206                for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
207                    j=pattern->mainPattern->index[iptr_ij];
208                    m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
209                    k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];
210                    ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];
211                    flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));
212                }
213                #pragma ivdep
214                for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
215                    j=pattern->couplePattern->index[iptr_ij];
216                    m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
217                    k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];
218                    ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];
219                    flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));
220                }
221             }
222          } else {
223             #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij)
224             for (i = 0; i < n; ++i) {
225                u_i=u[i];
226                #pragma ivdep
227                for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
228                    j=pattern->mainPattern->index[iptr_ij];
229                    m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
230                    flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*a*m_ij;
231                }
232                #pragma ivdep
233                for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
234                    j=pattern->couplePattern->index[iptr_ij];
235                    m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
236                    flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*a*m_ij;
237                }
238             }
239    
 void Paso_Solver_FluxControl_addDiffusion(Paso_Solver_FluxControl * fc, double alpha, Paso_SystemMatrix * B) {  
   dim_t n,i;  
   index_t color, iptr_ij,j,iptr_ji;  
   register double d_ij;  
240    
241    if (fc==NULL) return;        }
242    n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);    } else {
243    /* TODO test - same pattern + block size */        if ( (ABS(b) >0 ) ) {
244             #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,k_ij, ml_ij)
245             for (i = 0; i < n; ++i) {
246                u_i=u[i];
247                #pragma ivdep
248                for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
249                    j=pattern->mainPattern->index[iptr_ij];
250                    k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];
251                    ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];
252                    flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*b2*(ml_ij+k_ij);
253                }
254                #pragma ivdep
255                for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
256                    j=pattern->couplePattern->index[iptr_ij];
257                    k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];
258                    ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];
259                    flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*b2*(ml_ij+k_ij);
260                }
261             }
262    
263    #pragma omp parallel private(color)        }
264    {    }
265         /* process main block */  }
266         for (color=0;color<fc->num_colors;++color) {  /*
267             #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)   *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
268             for (i = 0; i < n; ++i) {   *  
269                 if (fc->colorOf[i]==color) {   */
270                    for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {  void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u)
271                       j=fc->matrix->mainBlock->pattern->index[iptr_ij];  {
272                       if (i<j) {    dim_t n, i, j;
273                          /* find entry a[j,i] */    Paso_SystemMatrixPattern *pattern;
274                          for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {    double *remote_u;
275                              if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {    register double u_i, f_ij;
276                                  d_ij=(-alpha)*MIN3(0.,fc->matrix->mainBlock->val[iptr_ij],fc->matrix->mainBlock->val[iptr_ji]);    register index_t iptr_ij;
277                                  B->mainBlock->val[iptr_ij]+=alpha*fc->matrix->mainBlock->val[iptr_ij]+d_ij;  
278                                  B->mainBlock->val[iptr_ji]+=alpha*fc->matrix->mainBlock->val[iptr_ji]+d_ij;    n=Paso_SystemMatrix_getTotalNumRows(f);
279                                  B->mainBlock->val[fc->main_iptr[i]]-=d_ij;    pattern=f->pattern;
280                                  B->mainBlock->val[fc->main_iptr[j]]-=d_ij;    remote_u=Paso_SystemMatrix_finishCollect(f);
281                                  break;    #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
282                              }    for (i = 0; i < n; ++i) {
283                          }        u_i=u[i];
284                       }        #pragma ivdep
285                                  for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
286                    }            j=pattern->mainPattern->index[iptr_ij];
287                    /* TODO process couple block */            f_ij=f->mainBlock->val[iptr_ij];
288                 }            if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
289             }        }
290             #pragma omp barrier        #pragma ivdep
291         }        for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
292              j=pattern->couplePattern->index[iptr_ij];
293              f_ij=f->coupleBlock->val[iptr_ij];
294              if (f_ij * (u_i-remote_u[j]) <= 0) f->coupleBlock->val[iptr_ij]=0;
295          }
296    }    }
   
297  }  }
 /**************************************************************/  
   
 /* adds antidiffusion to fa    
   
    P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])    
    P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])    
    Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])    
    Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])    
    d_ij=max(0,-a[i,j],-a[j,i])  
    l_ji=max(0,a[j,i],a[j,i]-a[i,j])  
    if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)  
       r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]  
       f_ij =min(FLUX_LIMITER(r_ij)*d_ij,l_ji) (u[i]-u[j])=min(FLUX_LIMITER(r_ij)*a[i,j],a[j,i]-a[i,j]) (u[i]-u[j])  
       fa[i]+=f_ij  
       fa[j]-=f_ij  
   
 */  
298    
 void Paso_Solver_FluxControl_setAntiDiffusiveFlux(Paso_Solver_FluxControl * fc, double * u, double* fa) {  
299    
300    register double u_i,P_p,P_n,Q_p,Q_n,r_p,r_n, a_ij, d, u_j, r_ij, f_ij, a_ji;  void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
301    double *u_remote=NULL;                                     const double* QN,const double* QP,double* RN,double* RP)
302    index_t color, iptr_ij,j,iptr_ji, i;  {
303    dim_t n;    dim_t n, i, j;
304      Paso_SystemMatrixPattern *pattern;
305      register double f_ij, PP_i, PN_i;
306      register index_t iptr_ij;
307    
308      n=Paso_SystemMatrix_getTotalNumRows(f);
309      pattern=f->pattern;
310      #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
311      for (i = 0; i < n; ++i) {
312          PP_i=0.;
313          PN_i=0.;
314          #pragma ivdep
315          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
316              j=pattern->mainPattern->index[iptr_ij];
317              if (i != j ) {
318                 f_ij=f->mainBlock->val[iptr_ij];
319                 if (f_ij <=0  ) {
320                    PN_i+=f_ij;
321                 } else {
322                    PP_i+=f_ij;
323                 }
324              }
325          }
326          #pragma ivdep
327          for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
328              j=pattern->couplePattern->index[iptr_ij];
329              f_ij=f->coupleBlock->val[iptr_ij];
330              if (f_ij <=0  ) {
331                    PN_i+=f_ij;
332              } else {
333                    PP_i+=f_ij;
334              }
335          }
336          if (PN_i<0) {
337             RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i);
338          } else {
339             RN[i]=1.;
340          }
341          if (PP_i>0) {
342             RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i);
343          } else {
344             RP[i]=1.;
345          }
346      }
347    
348    }
349    
350    if (fc==NULL) return;  void Paso_FCTransportProblem_addCorrectedFluxes(double* f,Paso_SystemMatrix *flux_matrix,const double* RN,const double* RP)
351    n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);  {
352    /* exchange */    dim_t n, i, j;
353    Paso_SystemMatrix_startCollect(fc->matrix,u);    Paso_SystemMatrixPattern *pattern;
354    u_remote=Paso_SystemMatrix_finishCollect(fc->matrix);    double *remote_RN, *remote_RP;
355      register double RN_i, RP_i, f_i, f_ij;
356    #pragma omp parallel private(color)    register index_t iptr_ij;
357    {    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
358         for (color=0;color<fc->num_colors;++color) {  
359             #pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,r_p,r_n,iptr_ij,a_ij,d,j,iptr_ji, u_j, r_ij, f_ij, a_ji)    Paso_SystemMatrix_startCollect(flux_matrix,RN);
360             for (i = 0; i < n; ++i) {    pattern=flux_matrix->pattern;
361                if (fc->colorOf[i]==color) {    remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix);
362                    u_i=u[i];    #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
363                    /* gather the smoothness sensor */    for (i = 0; i < n; ++i) {
364                    P_p=0.;       RN_i=RN[i];
365                    P_n=0.;       RP_i=RP[i];
366                    Q_p=0.;       f_i=0;
367                    Q_n=0.;       #pragma ivdep
368                    #pragma ivdep       for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
369                for (iptr_ij=(fc->matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {           j=pattern->mainPattern->index[iptr_ij];
370                        a_ij=fc->matrix->mainBlock->val[iptr_ij];           f_ij=flux_matrix->mainBlock->val[iptr_ij];
371                        j=fc->matrix->mainBlock->pattern->index[iptr_ij];           if (f_ij >=0) {
372                        d=u[j]-u_i;                f_i+=f_ij*MIN(RP_i,RN[j]);
373                        if (a_ij<0.) {           } else {
374                           if (d<0.) {                f_i+=f_ij*MIN(RN_i,RP[j]);
375                              P_p+=a_ij*d;           }
376                           } else {       }
377                              P_n+=a_ij*d;       #pragma ivdep
378                           }       for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
379                        } else {            j=pattern->couplePattern->index[iptr_ij];
380                           if (d>0.) {            f_ij=flux_matrix->coupleBlock->val[iptr_ij];
381                              Q_p+=a_ij*d;            if (f_ij >=0) {
382                           } else {                f_i+=f_ij*MIN(RP_i,remote_RN[j]);
383                              Q_n+=a_ij*d;            }
384                           }        }
385                        }        f[i]=f_i;
386                }    }
387                    #pragma ivdep    Paso_SystemMatrix_startCollect(flux_matrix,RP);
388                for (iptr_ij=(fc->matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {    remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);
389                        a_ij=fc->matrix->coupleBlock->val[iptr_ij];    #pragma omp parallel for schedule(static) private(i, RN_i, iptr_ij, j, f_ij, f_i)
390                        j=fc->matrix->coupleBlock->pattern->index[iptr_ij];    for (i = 0; i < n; ++i) {
391                        d=u_remote[j]-u_i;       RN_i=RN[i];
392                        if (a_ij<0.) {       f_i=0;
393                           if (d<0.) {       #pragma ivdep
394                              P_p+=a_ij*d;       for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
395                           } else {            j=pattern->couplePattern->index[iptr_ij];
396                              P_n+=a_ij*d;            f_ij=flux_matrix->coupleBlock->val[iptr_ij];
397                           }            if (f_ij < 0) {
398                        } else {                f_i+=f_ij*MIN(RN_i,remote_RP[j]);
399                           if (d>0.) {            }
400                              Q_p+=a_ij*d;       }
401                           } else {       f[i]+=f_i;
                             Q_n+=a_ij*d;  
                          }  
                       }  
               }  
                   /* set the smoothness indicators */  
                   r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;  
                   r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;  
                   /* anti diffusive flux from main block */  
                   for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {  
                      a_ij=fc->matrix->mainBlock->val[iptr_ij];  
                      j=fc->matrix->mainBlock->pattern->index[iptr_ij];  
                      if (a_ij < 0 && i!=j) {  
                         /* find entry a[j,i] */  
                         for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {  
                             if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {  
                                 a_ji=fc->matrix->mainBlock->val[iptr_ji];  
                                 if  (a_ji > a_ij || (a_ji == a_ij && j<i) ) {  
                                     u_j=u[j];  
                                     r_ij = u_i>u_j ? r_p : r_n;  
                                     f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j);  
                                     fa[i]+=f_ij;  
                                     fa[j]-=f_ij;  
                                     break;  
                                 }  
                             }  
                         }  
                      }  
                   }  
                   /* anti diffusive flux from couple block */  
   
                   /* TODO */  
               }  
            }  
            #pragma omp barrier  
        }  
402    }    }
403  }  }

Legend:
Removed from v.1361  
changed lines
  Added in v.1408

  ViewVC Help
Powered by ViewVC 1.1.26