/[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

revision 1370 by gross, Wed Jan 2 09:21:43 2008 UTC revision 1410 by gross, Thu Feb 7 04:24:00 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  /* adds A plus stabelising diffusion into the matrix B        */  /* 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  /* d_ij=alpha*max(0,-a[i,j],-a[j,i])  */   * if fc->iteration_matrix==NULL, fc->iteration_matrix is allocated
35  /* b[i,j]+=alpha*(a[i,j]+d_ij)  */   *
36  /* b[j,i]+=alpha*(a[j,i]+d_ij)  */   * a=transport_matrix
37  /* b[i,i]-=alpha*d_ij  */   * b= low_order_transport_matrix = - iteration_matrix
38  /* b[j,j]-=alpha*d_ij  */   * c=main diagonal low_order_transport_matrix
39     * initialize c[i] mit a[i,i]
40     *    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_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {  void Paso_FCTransportProblem_setLowOrderOperator(Paso_FCTransportProblem * fc) {
46    dim_t n,i;    dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i;
47    index_t color, iptr_ij,j,iptr_ji;    index_t color, iptr_ij,j,iptr_ji;
48    register double d_ij;    Paso_SystemMatrixPattern *pattern;
49      register double d_ij, sum, rtmp1, rtmp2;
50    
51    if (fc==NULL) return;    if (fc==NULL) return;
52    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);    if (fc->iteration_matrix==NULL) {
53          fc->iteration_matrix=Paso_SystemMatrix_alloc(fc->transport_matrix->type,
54    #pragma omp parallel private(color)                                                     fc->transport_matrix->pattern,
55    {                                                     fc->transport_matrix->row_block_size,
56         /* process main block */                                                     fc->transport_matrix->col_block_size);
        for (color=0;color<fc->num_colors;++color) {  
            #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)  
            for (i = 0; i < n; ++i) {  
                if (fc->colorOf[i]==color) {  
                   fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];  
   
                   for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {  
                      j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];  
                      if (j<i) {  
                         /* find entry a[j,i] */  
                         for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {  
                             if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {  
                                 d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);  
 printf("%d %d -> %e\n",i,j,d_ij);  
                                 fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;  
                                 fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij;  
 printf("%d %d -> %e -> %e %e \n",i,j,d_ij,fc->transport_matrix->mainBlock->val[iptr_ij],fc->transport_matrix->mainBlock->val[iptr_ji]);  
                                 fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;  
                                 fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;  
                                 break;  
                             }  
                         }  
                      }  
                           
                   }  
                   /* TODO process couple block */  
                }  
            }  
            #pragma omp barrier  
        }  
57    }    }
58    
59      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) {
64              sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
65              /* look at a[i,j] */
66              for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
67                  j=pattern->mainPattern->index[iptr_ij];
68                  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;
78                        }
79                     }
80                 }
81             }
82             /* TODO process couple block */
83            
84             /* set main diagonal entry */
85             fc->main_diagonal_low_order_transport_matrix[i]=sum;
86          }
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  void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {    n=Paso_SystemMatrix_getTotalNumRows(L);
108    
109    double *remote_u=NULL;    if (ABS(a)>0) {
110            Paso_SystemMatrix_startCollect(L,u);
111    if (fc==NULL) return;    }
112    Paso_SystemMatrix_startCollect(fc->transport_matrix,u);    #pragma omp parallel for private(i) schedule(static)
113    /* process main block */    for (i = 0; i < n; ++i) {
114    Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);        out[i]=M[i]*u[i]+b*Q[i];
115    /* finish exchange */    }
116    remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);    if (ABS(a)>0) {
117    /* process couple block */        pattern=L->pattern;
118    Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);        remote_u=Paso_SystemMatrix_finishCollect(L);
119          #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
120    Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);        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  /* adds antidiffusion to fa     *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
142     *
    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  
   
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=MAX(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=MAX(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  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {  /*
181     *
182    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;   *  f_{ij} = (m_{ij} - dt (1-theta) d_{ij}) (u_last[j]-u_last[i]) - (m_{ij} + dt theta d_{ij}) (u[j]-u[i])
183    index_t color, iptr_ij,j,iptr_ji, i;   *        
184    dim_t n;   * 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_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const double* u, const double* u_last)
188    {
189      dim_t n, j, i;
190      double *remote_u=NULL, *remote_u_last=NULL;
191      index_t iptr_ij;
192      const double f1=- dt * (1.-fc->theta);
193      const double f2=  dt * fc->theta;
194      register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;
195      Paso_SystemMatrixPattern *pattern;
196      Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);
197      Paso_SystemMatrix_startCollect(fc->iteration_matrix,u_last);
198      n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
199      pattern=fc->iteration_matrix->pattern;
200      remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
201      remote_u_last=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
202      if ( (ABS(f1) >0 ) ) {
203         if ( (ABS(f2) >0 ) ) {
204            #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
205            for (i = 0; i < n; ++i) {
206               u_i=u[i];
207               u_last_i=u_last[i];
208               #pragma ivdep
209               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
210                  j=pattern->mainPattern->index[iptr_ij];
211                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
212                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
213                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
214                  d_u=u[j]-u_i;
215                  d_u_last=u_last[j]-u_last_i;
216                  flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
217               }
218               #pragma ivdep
219               for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
220                  j=pattern->couplePattern->index[iptr_ij];
221                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
222                  d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
223                                                 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
224                  d_u=remote_u[j]-u_i;
225                  d_u_last=remote_u_last[j]-u_last_i;
226                  flux_matrix->coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
227               }
228            }
229         } else {
230            #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
231            for (i = 0; i < n; ++i) {
232               u_i=u[i];
233               u_last_i=u_last[i];
234               #pragma ivdep
235               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
236                  j=pattern->mainPattern->index[iptr_ij];
237                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
238                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
239                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
240                  d_u=u[j]-u_i;
241                  d_u_last=u_last[j]-u_last_i;
242                  flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
243               }
244               #pragma ivdep
245               for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
246                  j=pattern->couplePattern->index[iptr_ij];
247                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
248                  d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
249                                                 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
250                  d_u=remote_u[j]-u_i;
251                  d_u_last=remote_u_last[j]-u_last_i;
252                  flux_matrix->coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
253               }
254            }
255         }
256      } else {
257         if ( (ABS(f2) >0 ) ) {
258            #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
259            for (i = 0; i < n; ++i) {
260               u_i=u[i];
261               u_last_i=u_last[i];
262               #pragma ivdep
263               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
264                  j=pattern->mainPattern->index[iptr_ij];
265                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
266                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
267                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
268                  d_u=u[j]-u_i;
269                  d_u_last=u_last[j]-u_last_i;
270                  flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
271               }
272               #pragma ivdep
273               for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
274                  j=pattern->couplePattern->index[iptr_ij];
275                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
276                  d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
277                                                 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
278                  d_u=remote_u[j]-u_i;
279                  d_u_last=remote_u_last[j]-u_last_i;
280                  flux_matrix->coupleBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
281               }
282            }
283         } else {
284            #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
285            for (i = 0; i < n; ++i) {
286               u_i=u[i];
287               u_last_i=u_last[i];
288               #pragma ivdep
289               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
290                  j=pattern->mainPattern->index[iptr_ij];
291                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
292                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
293                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
294                  d_u=u[j]-u_i;
295                  d_u_last=u_last[j]-u_last_i;
296                  flux_matrix->mainBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
297               }
298               #pragma ivdep
299               for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
300                  j=pattern->couplePattern->index[iptr_ij];
301                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
302                  d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
303                                                 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
304                  d_u=remote_u[j]-u_i;
305                  d_u_last=remote_u_last[j]-u_last_i;
306                  flux_matrix->coupleBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
307               }
308            }
309          }
310      }
311    }
312    /*
313     *
314     * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])
315     *
316     * m=fc->mass matrix
317     * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
318     */
319    void Paso_FCTransportProblem_updateAntiDiffusionFlux(const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix,const double a, const double b, const double* u)
320    {
321      dim_t n, j, i;
322      double *remote_u;
323      index_t iptr_ij;
324      const double b2=-b;
325      register double m_ij, ml_ij, k_ij, u_i;
326      Paso_SystemMatrixPattern *pattern;
327      Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);
328      n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
329      pattern=fc->iteration_matrix->pattern;
330      remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
331      if ( (ABS(a) >0 ) ) {
332          if ( (ABS(b) >0 ) ) {
333             #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij)
334             for (i = 0; i < n; ++i) {
335                u_i=u[i];
336                #pragma ivdep
337                for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
338                    j=pattern->mainPattern->index[iptr_ij];
339                    m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
340                    k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];
341                    ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];
342                    flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));
343                }
344                #pragma ivdep
345                for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
346                    j=pattern->couplePattern->index[iptr_ij];
347                    m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
348                    k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];
349                    ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];
350                    flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));
351                }
352             }
353          } else {
354             #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij)
355             for (i = 0; i < n; ++i) {
356                u_i=u[i];
357                #pragma ivdep
358                for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
359                    j=pattern->mainPattern->index[iptr_ij];
360                    m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
361                    flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*a*m_ij;
362                }
363                #pragma ivdep
364                for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
365                    j=pattern->couplePattern->index[iptr_ij];
366                    m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
367                    flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*a*m_ij;
368                }
369             }
370    
371    
372          }
373      } else {
374          if ( (ABS(b) >0 ) ) {
375             #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,k_ij, ml_ij)
376             for (i = 0; i < n; ++i) {
377                u_i=u[i];
378                #pragma ivdep
379                for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
380                    j=pattern->mainPattern->index[iptr_ij];
381                    k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];
382                    ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];
383                    flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*b2*(ml_ij+k_ij);
384                }
385                #pragma ivdep
386                for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
387                    j=pattern->couplePattern->index[iptr_ij];
388                    k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];
389                    ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];
390                    flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*b2*(ml_ij+k_ij);
391                }
392             }
393    
394          }
395      }
396    }
397    /*
398     *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
399     *  
400     */
401    void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u)
402    {
403      dim_t n, i, j;
404      Paso_SystemMatrixPattern *pattern;
405      double *remote_u;
406      register double u_i, f_ij;
407      register index_t iptr_ij;
408    
409      n=Paso_SystemMatrix_getTotalNumRows(f);
410      pattern=f->pattern;
411      remote_u=Paso_SystemMatrix_finishCollect(f);
412      #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
413      for (i = 0; i < n; ++i) {
414          u_i=u[i];
415          #pragma ivdep
416          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
417              j=pattern->mainPattern->index[iptr_ij];
418              f_ij=f->mainBlock->val[iptr_ij];
419              if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
420          }
421          #pragma ivdep
422          for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
423              j=pattern->couplePattern->index[iptr_ij];
424              f_ij=f->coupleBlock->val[iptr_ij];
425              if (f_ij * (u_i-remote_u[j]) <= 0) f->coupleBlock->val[iptr_ij]=0;
426          }
427      }
428    }
429    
   if (fc==NULL) return;  
   n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);  
   /* exchange */  
430    
431    void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
432                                       const double* QN,const double* QP,double* RN,double* RP)
433    {
434      dim_t n, i, j;
435      Paso_SystemMatrixPattern *pattern;
436      register double f_ij, PP_i, PN_i;
437      register index_t iptr_ij;
438    
439      n=Paso_SystemMatrix_getTotalNumRows(f);
440      pattern=f->pattern;
441      #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
442      for (i = 0; i < n; ++i) {
443          PP_i=0.;
444          PN_i=0.;
445          #pragma ivdep
446          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
447              j=pattern->mainPattern->index[iptr_ij];
448              if (i != j ) {
449                 f_ij=f->mainBlock->val[iptr_ij];
450                 if (f_ij <=0) {
451                    PN_i+=f_ij;
452                 } else {
453                    PP_i+=f_ij;
454                 }
455              }
456          }
457          #pragma ivdep
458          for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
459              j=pattern->couplePattern->index[iptr_ij];
460              f_ij=f->coupleBlock->val[iptr_ij];
461              if (f_ij <=0  ) {
462                    PN_i+=f_ij;
463              } else {
464                    PP_i+=f_ij;
465              }
466          }
467          if (PN_i<0) {
468             RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i);
469          } else {
470             RN[i]=1.;
471          }
472          if (PP_i>0) {
473             RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i);
474          } else {
475             RP[i]=1.;
476          }
477      }
478    
479    #pragma omp parallel private(color)  }
   {  
        for (color=0;color<fc->num_colors;++color) {  
            #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)  
            for (i = 0; i < n; ++i) {  
               if (fc->colorOf[i]==color) {  
                   u_i=u[i];  
                   /* gather the smoothness sensor */  
                   P_p=0.;  
                   P_n=0.;  
                   Q_p=0.;  
                   Q_n=0.;  
                   #pragma ivdep  
               for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {  
                       a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];  
                       j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];  
                       d=u[j]-u_i;  
                       if (a_ij<0.) {  
                          if (d<0.) {  
                             P_p+=a_ij*d;  
                          } else {  
                             P_n+=a_ij*d;  
                          }  
                       } else {  
                          if (d>0.) {  
                             Q_p+=a_ij*d;  
                          } else {  
                             Q_n+=a_ij*d;  
                          }  
                       }  
               }  
                   #pragma ivdep  
               for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {  
                       a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];  
                       j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];  
                       d=u_remote[j]-u_i;  
                       if (a_ij<0.) {  
                          if (d<0.) {  
                             P_p+=a_ij*d;  
                          } else {  
                             P_n+=a_ij*d;  
                          }  
                       } else {  
                          if (d>0.) {  
                             Q_p+=a_ij*d;  
                          } else {  
                             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->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {  
                      a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];  
                      j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];  
                      if (a_ij < 0 && i!=j) {  
                         /* find entry a[j,i] */  
                         for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {  
                             if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {  
                                 a_ji=fc->flux_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 */  
480    
481                    /* TODO */  void Paso_FCTransportProblem_addCorrectedFluxes(double* f,Paso_SystemMatrix *flux_matrix,const double* RN,const double* RP)
482                }  {
483             }    dim_t n, i, j;
484             #pragma omp barrier    Paso_SystemMatrixPattern *pattern;
485         }    double *remote_RN, *remote_RP;
486      register double RN_i, RP_i, f_i, f_ij;
487      register index_t iptr_ij;
488      n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
489    
490      Paso_SystemMatrix_startCollect(flux_matrix,RN);
491      pattern=flux_matrix->pattern;
492      remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix);
493      #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
494      for (i = 0; i < n; ++i) {
495         RN_i=RN[i];
496         RP_i=RP[i];
497         f_i=0;
498         #pragma ivdep
499         for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
500             j=pattern->mainPattern->index[iptr_ij];
501             f_ij=flux_matrix->mainBlock->val[iptr_ij];
502             if (f_ij >=0) {
503                  f_i+=f_ij*MIN(RP_i,RN[j]);
504             } else {
505                  f_i+=f_ij*MIN(RN_i,RP[j]);
506             }
507         }
508         #pragma ivdep
509         for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
510              j=pattern->couplePattern->index[iptr_ij];
511              f_ij=flux_matrix->coupleBlock->val[iptr_ij];
512              if (f_ij >=0) {
513                  f_i+=f_ij*MIN(RP_i,remote_RN[j]);
514              }
515          }
516          f[i]+=f_i;
517      }
518      Paso_SystemMatrix_startCollect(flux_matrix,RP);
519      remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);
520      #pragma omp parallel for schedule(static) private(i, RN_i, iptr_ij, j, f_ij, f_i)
521      for (i = 0; i < n; ++i) {
522         RN_i=RN[i];
523         f_i=0;
524         #pragma ivdep
525         for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
526              j=pattern->couplePattern->index[iptr_ij];
527              f_ij=flux_matrix->coupleBlock->val[iptr_ij];
528              if (f_ij < 0) {
529                  f_i+=f_ij*MIN(RN_i,remote_RP[j]);
530              }
531         }
532         f[i]+=f_i;
533    }    }
534  }  }

Legend:
Removed from v.1370  
changed lines
  Added in v.1410

  ViewVC Help
Powered by ViewVC 1.1.26