/[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 1408 by gross, Mon Feb 4 07:19:50 2008 UTC revision 1562 by gross, Wed May 21 13:04:40 2008 UTC
# Line 63  void Paso_FCTransportProblem_setLowOrder Line 63  void Paso_FCTransportProblem_setLowOrder
63        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
64            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
65            /* look at a[i,j] */            /* look at a[i,j] */
66              #pragma ivdep
67            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
68                j=pattern->mainPattern->index[iptr_ij];                j=pattern->mainPattern->index[iptr_ij];
69                if (j!=i) {                if (j!=i) {
# Line 79  void Paso_FCTransportProblem_setLowOrder Line 80  void Paso_FCTransportProblem_setLowOrder
80                   }                   }
81               }               }
82           }           }
83           /* TODO process couple block */            #pragma ivdep
84                      for (iptr_ij=pattern->col_couplePattern->ptr[i];iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
85                  j=pattern->col_couplePattern->index[iptr_ij];
86                  /* find entry a[j,i] */
87                  for (iptr_ji=pattern->row_couplePattern->ptr[j]; iptr_ji<pattern->row_couplePattern->ptr[j+1]; ++iptr_ji) {
88                        if (pattern->row_couplePattern->index[iptr_ji]==i) {
89                            rtmp1=fc->transport_matrix->col_coupleBlock->val[iptr_ij];
90                            rtmp2=fc->transport_matrix->row_coupleBlock->val[iptr_ji];
91                            d_ij=-MIN3(0.,rtmp1,rtmp2);
92                            fc->iteration_matrix->col_coupleBlock->val[iptr_ij]=-(rtmp1+d_ij);
93                            fc->iteration_matrix->row_coupleBlock->val[iptr_ji]=-(rtmp2+d_ij);
94                            sum-=d_ij;
95                            break;
96                        }
97                  }
98             }
99           /* set main diagonal entry */           /* set main diagonal entry */
100           fc->main_diagonal_low_order_transport_matrix[i]=sum;           fc->main_diagonal_low_order_transport_matrix[i]=sum;
101        }        }
# Line 92  void Paso_FCTransportProblem_setLowOrder Line 107  void Paso_FCTransportProblem_setLowOrder
107   */   */
108  void Paso_SolverFCT_setMuPaLuPbQ(double* out,  void Paso_SolverFCT_setMuPaLuPbQ(double* out,
109                                   const double* M,                                   const double* M,
110                                   const double* u,                                   const Paso_Coupler* u_coupler,
111                                   const double a,                                   const double a,
112                                   Paso_SystemMatrix *L,                                   const Paso_SystemMatrix *L,
113                                   const double b,                                   const double b,
114                                   const double* Q)                                   const double* Q)
115  {  {
116    dim_t n, i, j;    dim_t n, i, j;
117    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
118    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
119      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
120    register double sum, u_i, l_ij;    register double sum, u_i, l_ij;
121    register index_t iptr_ij;    register index_t iptr_ij;
122    
123    n=Paso_SystemMatrix_getTotalNumRows(L);    n=Paso_SystemMatrix_getTotalNumRows(L);
124    
   if (ABS(a)>0) {  
       Paso_SystemMatrix_startCollect(L,u);  
   }  
125    #pragma omp parallel for private(i) schedule(static)    #pragma omp parallel for private(i) schedule(static)
126    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
127        out[i]=M[i]*u[i]+b*Q[i];        out[i]=M[i]*u[i]+b*Q[i];
128    }    }
129    if (ABS(a)>0) {    if (ABS(a)>0) {
130        pattern=L->pattern;        pattern=L->pattern;
       remote_u=Paso_SystemMatrix_finishCollect(L);  
131        #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)        #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
132        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
133            sum=0;            sum=0;
# Line 127  void Paso_SolverFCT_setMuPaLuPbQ(double* Line 139  void Paso_SolverFCT_setMuPaLuPbQ(double*
139                 sum+=l_ij*(u[j]-u_i);                 sum+=l_ij*(u[j]-u_i);
140            }            }
141            #pragma ivdep            #pragma ivdep
142        for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {        for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
143                 j=pattern->couplePattern->index[iptr_ij];                 j=pattern->col_couplePattern->index[iptr_ij];
144                 l_ij=L->coupleBlock->val[iptr_ij];                 l_ij=L->col_coupleBlock->val[iptr_ij];
145                 sum+=l_ij*(remote_u[j]-u_i);                 sum+=l_ij*(remote_u[j]-u_i);
146            }            }
147            out[i]+=a*sum;            out[i]+=a*sum;
# Line 141  void Paso_SolverFCT_setMuPaLuPbQ(double* Line 153  void Paso_SolverFCT_setMuPaLuPbQ(double*
153   *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )   *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
154   *   *
155  */  */
156  void Paso_SolverFCT_setQs(const double* u,double* QN, double* QP, Paso_SystemMatrix *L)  void Paso_SolverFCT_setQs(const Paso_Coupler* u_coupler,double* QN, double* QP, const Paso_SystemMatrix *L)
157  {  {
158    dim_t n, i, j;    dim_t n, i, j;
159    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
160    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
161      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
162    register double u_i, u_min_i, u_max_i, u_j;    register double u_i, u_min_i, u_max_i, u_j;
163    register index_t iptr_ij;    register index_t iptr_ij;
   
   Paso_SystemMatrix_startCollect(L,u);  
164    n=Paso_SystemMatrix_getTotalNumRows(L);    n=Paso_SystemMatrix_getTotalNumRows(L);
165    pattern=L->pattern;    pattern=L->pattern;
   remote_u=Paso_SystemMatrix_finishCollect(L);  
166    #pragma omp parallel for schedule(static) private(i, u_i, u_min_i, u_max_i, j, u_j, iptr_ij)    #pragma omp parallel for schedule(static) private(i, u_i, u_min_i, u_max_i, j, u_j, iptr_ij)
167    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
168       u_i=u[i];       u_i=u[i];
# Line 163  void Paso_SolverFCT_setQs(const double* Line 173  void Paso_SolverFCT_setQs(const double*
173           j=pattern->mainPattern->index[iptr_ij];           j=pattern->mainPattern->index[iptr_ij];
174           u_j=u[j];           u_j=u[j];
175           u_min_i=MIN(u_min_i,u_j);           u_min_i=MIN(u_min_i,u_j);
176           u_max_i=MIN(u_max_i,u_j);           u_max_i=MAX(u_max_i,u_j);
177       }       }
178       #pragma ivdep       #pragma ivdep
179       for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {       for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
180            j=pattern->couplePattern->index[iptr_ij];            j=pattern->col_couplePattern->index[iptr_ij];
181            u_j=remote_u[j];            u_j=remote_u[j];
182            u_min_i=MIN(u_min_i,u_j);            u_min_i=MIN(u_min_i,u_j);
183            u_max_i=MIN(u_max_i,u_j);            u_max_i=MAX(u_max_i,u_j);
184        }        }
185        QN[i]=u_min_i-u_i;        QN[i]=u_min_i-u_i;
186        QP[i]=u_max_i-u_i;        QP[i]=u_max_i-u_i;
# Line 179  void Paso_SolverFCT_setQs(const double* Line 189  void Paso_SolverFCT_setQs(const double*
189    
190  /*  /*
191   *   *
192   * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])   *  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])
193   *   *        
194   * m=fc->mass matrix   * m=fc->mass matrix
195   * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)   * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
196   */   */
197  void Paso_FCTransportProblem_updateAntiDiffusionFlux(const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix,const double a, const double b, const double* u)  void Paso_FCTransportProblem_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const Paso_Coupler* u_coupler)
198  {  {
199    dim_t n, j, i;    dim_t n, j, i;
   double *remote_u;  
200    index_t iptr_ij;    index_t iptr_ij;
201    const double b2=-b;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
202    register double m_ij, ml_ij, k_ij, u_i;    const double *u_last= Paso_Coupler_borrowLocalData(fc->u_coupler);
203    Paso_SystemMatrixPattern *pattern;    const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
204    Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);    const double *remote_u_last=Paso_Coupler_borrowRemoteData(fc->u_coupler);
205      const double f1=- dt * (1.-fc->theta);
206      const double f2=  dt * fc->theta;
207      register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;
208      const Paso_SystemMatrixPattern *pattern=fc->iteration_matrix->pattern;
209    n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
210    pattern=fc->iteration_matrix->pattern;    if ( (ABS(f1) >0 ) ) {
211    remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);       if ( (ABS(f2) >0 ) ) {
212            #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)
213    if ( (ABS(a) >0 ) ) {          for (i = 0; i < n; ++i) {
214        if ( (ABS(b) >0 ) ) {             u_i=u[i];
215           #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij)             u_last_i=u_last[i];
216           for (i = 0; i < n; ++i) {             #pragma ivdep
217              u_i=u[i];             for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
218              #pragma ivdep                j=pattern->mainPattern->index[iptr_ij];
219              for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {                m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
220                  j=pattern->mainPattern->index[iptr_ij];                d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
221                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];                                                  fc->iteration_matrix->mainBlock->val[iptr_ij]);
222                  k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];                d_u=u[j]-u_i;
223                  ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];                d_u_last=u_last[j]-u_last_i;
224                  flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));                flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
225              }             }
226              #pragma ivdep             #pragma ivdep
227              for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {             for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
228                  j=pattern->couplePattern->index[iptr_ij];                j=pattern->col_couplePattern->index[iptr_ij];
229                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];                m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
230                  k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];                d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
231                  ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];                                               fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
232                  flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));                d_u=remote_u[j]-u_i;
233              }                d_u_last=remote_u_last[j]-u_last_i;
234           }                flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
235        } else {             }
236           #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij)          }
237           for (i = 0; i < n; ++i) {       } else {
238              u_i=u[i];          #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)
239              #pragma ivdep          for (i = 0; i < n; ++i) {
240              for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {             u_i=u[i];
241                  j=pattern->mainPattern->index[iptr_ij];             u_last_i=u_last[i];
242                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];             #pragma ivdep
243                  flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*a*m_ij;             for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
244              }                j=pattern->mainPattern->index[iptr_ij];
245              #pragma ivdep                m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
246              for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {                d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
247                  j=pattern->couplePattern->index[iptr_ij];                                                  fc->iteration_matrix->mainBlock->val[iptr_ij]);
248                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];                d_u=u[j]-u_i;
249                  flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*a*m_ij;                d_u_last=u_last[j]-u_last_i;
250              }                flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
251           }             }
252               #pragma ivdep
253               for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
254        }                j=pattern->col_couplePattern->index[iptr_ij];
255                  m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
256                  d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
257                                                 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
258                  d_u=remote_u[j]-u_i;
259                  d_u_last=remote_u_last[j]-u_last_i;
260                  flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
261               }
262            }
263         }
264    } else {    } else {
265        if ( (ABS(b) >0 ) ) {       if ( (ABS(f2) >0 ) ) {
266           #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,k_ij, ml_ij)          #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)
267           for (i = 0; i < n; ++i) {          for (i = 0; i < n; ++i) {
268              u_i=u[i];             u_i=u[i];
269              #pragma ivdep             u_last_i=u_last[i];
270              for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {             #pragma ivdep
271                  j=pattern->mainPattern->index[iptr_ij];             for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
272                  k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];                j=pattern->mainPattern->index[iptr_ij];
273                  ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];                m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
274                  flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*b2*(ml_ij+k_ij);                d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
275              }                                                  fc->iteration_matrix->mainBlock->val[iptr_ij]);
276              #pragma ivdep                d_u=u[j]-u_i;
277              for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {                d_u_last=u_last[j]-u_last_i;
278                  j=pattern->couplePattern->index[iptr_ij];                flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
279                  k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];             }
280                  ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];             #pragma ivdep
281                  flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*b2*(ml_ij+k_ij);             for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
282              }                j=pattern->col_couplePattern->index[iptr_ij];
283           }                m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
284                  d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
285        }                                               fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
286                  d_u=remote_u[j]-u_i;
287                  d_u_last=remote_u_last[j]-u_last_i;
288                  flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
289               }
290            }
291         } else {
292            #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)
293            for (i = 0; i < n; ++i) {
294               u_i=u[i];
295               u_last_i=u_last[i];
296               #pragma ivdep
297               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
298                  j=pattern->mainPattern->index[iptr_ij];
299                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
300                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
301                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
302                  d_u=u[j]-u_i;
303                  d_u_last=u_last[j]-u_last_i;
304                  flux_matrix->mainBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
305               }
306               #pragma ivdep
307               for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
308                  j=pattern->col_couplePattern->index[iptr_ij];
309                  m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
310                  d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
311                                                 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
312                  d_u=remote_u[j]-u_i;
313                  d_u_last=remote_u_last[j]-u_last_i;
314                  flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
315               }
316            }
317          }
318    }    }
319  }  }
320  /*  /*
321   *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0   *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
322   *     *  
323   */   */
324  void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u)  void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const Paso_Coupler* u_coupler)
325  {  {
326    dim_t n, i, j;    dim_t n, i, j;
327    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
328    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
329      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
330    register double u_i, f_ij;    register double u_i, f_ij;
331    register index_t iptr_ij;    register index_t iptr_ij;
332    
333    n=Paso_SystemMatrix_getTotalNumRows(f);    n=Paso_SystemMatrix_getTotalNumRows(f);
334    pattern=f->pattern;    pattern=f->pattern;
   remote_u=Paso_SystemMatrix_finishCollect(f);  
335    #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)    #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
336    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
337        u_i=u[i];        u_i=u[i];
# Line 288  void Paso_FCTransportProblem_applyPreAnt Line 342  void Paso_FCTransportProblem_applyPreAnt
342            if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;            if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
343        }        }
344        #pragma ivdep        #pragma ivdep
345        for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {        for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
346            j=pattern->couplePattern->index[iptr_ij];            j=pattern->col_couplePattern->index[iptr_ij];
347            f_ij=f->coupleBlock->val[iptr_ij];            f_ij=f->col_coupleBlock->val[iptr_ij];
348            if (f_ij * (u_i-remote_u[j]) <= 0) f->coupleBlock->val[iptr_ij]=0;            if (f_ij * (u_i-remote_u[j]) <= 0) f->col_coupleBlock->val[iptr_ij]=0;
349        }        }
350    }    }
351  }  }
352    
353    
354  void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,  void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
355                                     const double* QN,const double* QP,double* RN,double* RP)                                     const Paso_Coupler* QN_coupler,const Paso_Coupler* QP_coupler,double* RN,double* RP)
356  {  {
357    dim_t n, i, j;    dim_t n, i, j;
358    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
359      const double *QN=Paso_Coupler_borrowLocalData(QN_coupler);
360      const double *QP=Paso_Coupler_borrowLocalData(QP_coupler);
361    register double f_ij, PP_i, PN_i;    register double f_ij, PP_i, PN_i;
362    register index_t iptr_ij;    register index_t iptr_ij;
363    
# Line 316  void Paso_FCTransportProblem_setRs(const Line 372  void Paso_FCTransportProblem_setRs(const
372            j=pattern->mainPattern->index[iptr_ij];            j=pattern->mainPattern->index[iptr_ij];
373            if (i != j ) {            if (i != j ) {
374               f_ij=f->mainBlock->val[iptr_ij];               f_ij=f->mainBlock->val[iptr_ij];
375               if (f_ij <=0  ) {               if (f_ij <=0) {
376                  PN_i+=f_ij;                  PN_i+=f_ij;
377               } else {               } else {
378                  PP_i+=f_ij;                  PP_i+=f_ij;
# Line 324  void Paso_FCTransportProblem_setRs(const Line 380  void Paso_FCTransportProblem_setRs(const
380            }            }
381        }        }
382        #pragma ivdep        #pragma ivdep
383        for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {        for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
384            j=pattern->couplePattern->index[iptr_ij];            j=pattern->col_couplePattern->index[iptr_ij];
385            f_ij=f->coupleBlock->val[iptr_ij];            f_ij=f->col_coupleBlock->val[iptr_ij];
386            if (f_ij <=0  ) {            if (f_ij <=0  ) {
387                  PN_i+=f_ij;                  PN_i+=f_ij;
388            } else {            } else {
# Line 347  void Paso_FCTransportProblem_setRs(const Line 403  void Paso_FCTransportProblem_setRs(const
403    
404  }  }
405    
406  void Paso_FCTransportProblem_addCorrectedFluxes(double* f,Paso_SystemMatrix *flux_matrix,const double* RN,const double* RP)  void Paso_FCTransportProblem_addCorrectedFluxes(double* f,const Paso_SystemMatrix *flux_matrix, const Paso_Coupler* RN_coupler, const Paso_Coupler* RP_coupler)
407  {  {
408    dim_t n, i, j;    dim_t n, i, j;
409    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
   double *remote_RN, *remote_RP;  
410    register double RN_i, RP_i, f_i, f_ij;    register double RN_i, RP_i, f_i, f_ij;
411    register index_t iptr_ij;    register index_t iptr_ij;
412      const double *RN=Paso_Coupler_borrowLocalData(RN_coupler);
413      const double *remote_RN=Paso_Coupler_borrowRemoteData(RN_coupler);
414      const double *RP=Paso_Coupler_borrowLocalData(RP_coupler);
415      const double *remote_RP=Paso_Coupler_borrowRemoteData(RP_coupler);
416    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
417    
   Paso_SystemMatrix_startCollect(flux_matrix,RN);  
418    pattern=flux_matrix->pattern;    pattern=flux_matrix->pattern;
   remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix);  
419    #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)    #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
420    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
421       RN_i=RN[i];       RN_i=RN[i];
# Line 375  void Paso_FCTransportProblem_addCorrecte Line 432  void Paso_FCTransportProblem_addCorrecte
432           }           }
433       }       }
434       #pragma ivdep       #pragma ivdep
435       for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {       for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
436            j=pattern->couplePattern->index[iptr_ij];            j=pattern->col_couplePattern->index[iptr_ij];
437            f_ij=flux_matrix->coupleBlock->val[iptr_ij];            f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];
438            if (f_ij >=0) {            if (f_ij >=0) {
439                f_i+=f_ij*MIN(RP_i,remote_RN[j]);                f_i+=f_ij*MIN(RP_i,remote_RN[j]);
440            }            }else {
       }  
       f[i]=f_i;  
   }  
   Paso_SystemMatrix_startCollect(flux_matrix,RP);  
   remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);  
   #pragma omp parallel for schedule(static) private(i, RN_i, iptr_ij, j, f_ij, f_i)  
   for (i = 0; i < n; ++i) {  
      RN_i=RN[i];  
      f_i=0;  
      #pragma ivdep  
      for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {  
           j=pattern->couplePattern->index[iptr_ij];  
           f_ij=flux_matrix->coupleBlock->val[iptr_ij];  
           if (f_ij < 0) {  
441                f_i+=f_ij*MIN(RN_i,remote_RP[j]);                f_i+=f_ij*MIN(RN_i,remote_RP[j]);
442            }            }
443       }        }
444       f[i]+=f_i;        f[i]+=f_i;
445    }    }
446  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26