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

revision 1401 by gross, Fri Jan 25 04:31:18 2008 UTC 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 MC(a,b) (FLUX_S(a,b)*MIN3(ABS((a)+(b))/2,2*ABS(a),2*ABS(b)))

#define FLUX_L(a,b) MC(a,b)  /* alter for other flux limiter */

#define FLUX_LIMITER(a) FLUX_L(1,a)
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, sum;    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                                                       fc->transport_matrix->pattern,
55                                                       fc->transport_matrix->row_block_size,
56                                                       fc->transport_matrix->col_block_size);
57      }
58
59    #pragma omp parallel for private(i,iptr_ij,j,iptr_ji,d_ij,sum)  schedule(static)    if (Paso_noError()) {
60    for (i = 0; i < n; ++i) {        pattern=fc->iteration_matrix->pattern;
61       sum=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];        n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
62       for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {        #pragma omp parallel for private(i,sum,iptr_ij,j,iptr_ji,rtmp1, rtmp2,d_ij)  schedule(static)
63           j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];        for (i = 0; i < n; ++i) {
64           if (j!=i) {            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
65               /* find entry a[j,i] */            /* look at a[i,j] */
66               for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
67                  if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {                j=pattern->mainPattern->index[iptr_ij];
68                      d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],                if (j!=i) {
69                                            fc->flux_matrix->mainBlock->val[iptr_ji]);                   /* find entry a[j,i] */
70                      fc->transport_matrix->mainBlock->val[iptr_ij]+=                   for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
71                                           alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;                      if (pattern->mainPattern->index[iptr_ji]==i) {
72                      sum-=d_ij;                          rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];
73                      break;                          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       /* TODO process couple block */
84             /* set main diagonal entry */
85       /* update main diagonal */           fc->main_diagonal_low_order_transport_matrix[i]=sum;
86       fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=sum;        }
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);
/*
*   sets fa=transport_matrix*u+anti-diffuison_flux(u)
*/

double *remote_u=NULL;

if (fc==NULL) return;
Paso_SystemMatrix_startCollect(fc->transport_matrix,u);
/* process main block */
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);
/* finish exchange */
remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);
/* process couple block */
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);
108
109    Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);    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  /* 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=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  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {  /*
181     *
182     * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])
183     *
184     * 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
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, d_ij, sum;
index_t color, iptr_ij,j,iptr_ji, i;
dim_t n;
240
241          }
242      } else {
243          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    if (fc==NULL) return;        }
264    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);    }
265    }
266    /*
267     *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
268     *
269     */
270    void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u)
271    {
272      dim_t n, i, j;
273      Paso_SystemMatrixPattern *pattern;
274      double *remote_u;
275      register double u_i, f_ij;
276      register index_t iptr_ij;
277
278      n=Paso_SystemMatrix_getTotalNumRows(f);
279      pattern=f->pattern;
280      remote_u=Paso_SystemMatrix_finishCollect(f);
281      #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              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 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    }
298
299
300    #pragma omp parallel  void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
301    {                                     const double* QN,const double* QP,double* RN,double* RP)
302        /*  {
303         * calculate the smootness sensors    dim_t n, i, j;
304        */    Paso_SystemMatrixPattern *pattern;
305        #pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,iptr_ij,a_ij,j,d)    register double f_ij, PP_i, PN_i;
306        for (i = 0; i < n; ++i) {    register index_t iptr_ij;
307            u_i=u[i];
308            P_p=0.;    n=Paso_SystemMatrix_getTotalNumRows(f);
309            P_n=0.;    pattern=f->pattern;
310            Q_p=0.;    #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
311            Q_n=0.;    for (i = 0; i < n; ++i) {
312            #pragma ivdep        PP_i=0.;
313        for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {        PN_i=0.;
314                 a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];        #pragma ivdep
315                 j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];        for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
316                 d=u[j]-u_i;            j=pattern->mainPattern->index[iptr_ij];
317                 if (a_ij<0.) {            if (i != j ) {
318                    if (d<0.) {               f_ij=f->mainBlock->val[iptr_ij];
319                        P_p+=a_ij*d;               if (f_ij <=0  ) {
320                    } else {                  PN_i+=f_ij;
321                        P_n+=a_ij*d;               } else {
322                    }                  PP_i+=f_ij;
323                 } 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 */
fc->r_p[i] = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
fc->r_n[i] = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.;

} /* end of row loop for smootheness indicator */

/*
* calculate antidiffusion
*/
#pragma omp for schedule(static) private(i, u_i, sum, iptr_ij, a_ij, j, iptr_ji, a_ji,d_ij, u_j, r_ij, f_ij)
for (i = 0; i < n; ++i) {
u_i=u[i];
sum=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 ( i!=j ) {
/* 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) {
a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
d_ij=-MIN3(0,a_ji,a_ij);
if (d_ij > 0.) {
u_j=u[j];
if (a_ji >= a_ij) {
r_ij = u_i>u_j ? fc->r_p[i] : fc->r_n[i];
f_ij =MIN(r_ij*d_ij,a_ji+d_ij);
} else {
r_ij = u_j>u_i ? fc->r_p[j] : fc->r_n[j];
f_ij =MIN(r_ij*d_ij,a_ij+d_ij);
}
sum+=f_ij*(u_i-u_j);
}
break;
}
}
}
324            }            }
325            /* anti diffusive flux from couple block */        }
326            /* TODO */        #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            fa[i]+=sum;  void Paso_FCTransportProblem_addCorrectedFluxes(double* f,Paso_SystemMatrix *flux_matrix,const double* RN,const double* RP)
351    {
352      dim_t n, i, j;
353      Paso_SystemMatrixPattern *pattern;
354      double *remote_RN, *remote_RP;
355      register double RN_i, RP_i, f_i, f_ij;
356      register index_t iptr_ij;
357      n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
358
359      Paso_SystemMatrix_startCollect(flux_matrix,RN);
360      pattern=flux_matrix->pattern;
361      remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix);
362      #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
363      for (i = 0; i < n; ++i) {
364         RN_i=RN[i];
365         RP_i=RP[i];
366         f_i=0;
367         #pragma ivdep
368         for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
369             j=pattern->mainPattern->index[iptr_ij];
370             f_ij=flux_matrix->mainBlock->val[iptr_ij];
371             if (f_ij >=0) {
372                  f_i+=f_ij*MIN(RP_i,RN[j]);
373             } else {
374                  f_i+=f_ij*MIN(RN_i,RP[j]);
375             }
376         }
377         #pragma ivdep
378         for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
379              j=pattern->couplePattern->index[iptr_ij];
380              f_ij=flux_matrix->coupleBlock->val[iptr_ij];
381              if (f_ij >=0) {
382                  f_i+=f_ij*MIN(RP_i,remote_RN[j]);
383              }
384        }        }
385    } /* end of parallel block */        f[i]=f_i;
386      }
387      Paso_SystemMatrix_startCollect(flux_matrix,RP);
388      remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);
389      #pragma omp parallel for schedule(static) private(i, RN_i, iptr_ij, j, f_ij, f_i)
390      for (i = 0; i < n; ++i) {
391         RN_i=RN[i];
392         f_i=0;
393         #pragma ivdep
394         for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
395              j=pattern->couplePattern->index[iptr_ij];
396              f_ij=flux_matrix->coupleBlock->val[iptr_ij];
397              if (f_ij < 0) {
398                  f_i+=f_ij*MIN(RN_i,remote_RP[j]);
399              }
400         }
401         f[i]+=f_i;
402      }
403  }  }

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