/[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 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1 
 /* $Id:$ */  
1    
2  /*******************************************************  /*******************************************************
3   *  *
4   *       Copyright 2007 by University of Queensland  * Copyright (c) 2003-2008 by University of Queensland
5   *  * Earth Systems Science Computational Center (ESSCC)
6   *                http://esscc.uq.edu.au  * http://www.uq.edu.au/esscc
7   *        Primary Business: Queensland, Australia  *
8   *  Licensed under the Open Software License version 3.0  * Primary Business: Queensland, Australia
9   *     http://www.opensource.org/licenses/osl-3.0.php  * Licensed under the Open Software License version 3.0
10   *  * http://www.opensource.org/licenses/osl-3.0.php
11   *******************************************************/  *
12    *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 44  Line 45 
45    
46  void Paso_FCTransportProblem_setLowOrderOperator(Paso_FCTransportProblem * fc) {  void Paso_FCTransportProblem_setLowOrderOperator(Paso_FCTransportProblem * fc) {
47    dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i;    dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i;
48    index_t color, iptr_ij,j,iptr_ji;    index_t iptr_ij,j,iptr_ji;
49    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
50    register double d_ij, sum, rtmp1, rtmp2;    register double d_ij, sum, rtmp1, rtmp2;
51    
# Line 63  void Paso_FCTransportProblem_setLowOrder Line 64  void Paso_FCTransportProblem_setLowOrder
64        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
65            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
66            /* look at a[i,j] */            /* look at a[i,j] */
67              #pragma ivdep
68            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) {
69                j=pattern->mainPattern->index[iptr_ij];                j=pattern->mainPattern->index[iptr_ij];
70                if (j!=i) {                if (j!=i) {
# Line 79  void Paso_FCTransportProblem_setLowOrder Line 81  void Paso_FCTransportProblem_setLowOrder
81                   }                   }
82               }               }
83           }           }
84           /* TODO process couple block */            #pragma ivdep
85                      for (iptr_ij=pattern->col_couplePattern->ptr[i];iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
86                  j=pattern->col_couplePattern->index[iptr_ij];
87                  /* find entry a[j,i] */
88                  for (iptr_ji=pattern->row_couplePattern->ptr[j]; iptr_ji<pattern->row_couplePattern->ptr[j+1]; ++iptr_ji) {
89                        if (pattern->row_couplePattern->index[iptr_ji]==i) {
90                            rtmp1=fc->transport_matrix->col_coupleBlock->val[iptr_ij];
91                            rtmp2=fc->transport_matrix->row_coupleBlock->val[iptr_ji];
92                            d_ij=-MIN3(0.,rtmp1,rtmp2);
93                            fc->iteration_matrix->col_coupleBlock->val[iptr_ij]=-(rtmp1+d_ij);
94                            fc->iteration_matrix->row_coupleBlock->val[iptr_ji]=-(rtmp2+d_ij);
95                            sum-=d_ij;
96                            break;
97                        }
98                  }
99             }
100           /* set main diagonal entry */           /* set main diagonal entry */
101           fc->main_diagonal_low_order_transport_matrix[i]=sum;           fc->main_diagonal_low_order_transport_matrix[i]=sum;
102        }        }
# Line 92  void Paso_FCTransportProblem_setLowOrder Line 108  void Paso_FCTransportProblem_setLowOrder
108   */   */
109  void Paso_SolverFCT_setMuPaLuPbQ(double* out,  void Paso_SolverFCT_setMuPaLuPbQ(double* out,
110                                   const double* M,                                   const double* M,
111                                   const double* u,                                   const Paso_Coupler* u_coupler,
112                                   const double a,                                   const double a,
113                                   Paso_SystemMatrix *L,                                   const Paso_SystemMatrix *L,
114                                   const double b,                                   const double b,
115                                   const double* Q)                                   const double* Q)
116  {  {
117    dim_t n, i, j;    dim_t n, i, j;
118    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
119    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
120      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
121    register double sum, u_i, l_ij;    register double sum, u_i, l_ij;
122    register index_t iptr_ij;    register index_t iptr_ij;
123    
124    n=Paso_SystemMatrix_getTotalNumRows(L);    n=Paso_SystemMatrix_getTotalNumRows(L);
125    
   if (ABS(a)>0) {  
       Paso_SystemMatrix_startCollect(L,u);  
   }  
126    #pragma omp parallel for private(i) schedule(static)    #pragma omp parallel for private(i) schedule(static)
127    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
128        out[i]=M[i]*u[i]+b*Q[i];        out[i]=M[i]*u[i]+b*Q[i];
129    }    }
130    if (ABS(a)>0) {    if (ABS(a)>0) {
131        pattern=L->pattern;        pattern=L->pattern;
       remote_u=Paso_SystemMatrix_finishCollect(L);  
132        #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)
133        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
134            sum=0;            sum=0;
# Line 127  void Paso_SolverFCT_setMuPaLuPbQ(double* Line 140  void Paso_SolverFCT_setMuPaLuPbQ(double*
140                 sum+=l_ij*(u[j]-u_i);                 sum+=l_ij*(u[j]-u_i);
141            }            }
142            #pragma ivdep            #pragma ivdep
143        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) {
144                 j=pattern->couplePattern->index[iptr_ij];                 j=pattern->col_couplePattern->index[iptr_ij];
145                 l_ij=L->coupleBlock->val[iptr_ij];                 l_ij=L->col_coupleBlock->val[iptr_ij];
146                 sum+=l_ij*(remote_u[j]-u_i);                 sum+=l_ij*(remote_u[j]-u_i);
147            }            }
148            out[i]+=a*sum;            out[i]+=a*sum;
# Line 141  void Paso_SolverFCT_setMuPaLuPbQ(double* Line 154  void Paso_SolverFCT_setMuPaLuPbQ(double*
154   *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )   *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
155   *   *
156  */  */
157  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)
158  {  {
159    dim_t n, i, j;    dim_t n, i, j;
160    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
161    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
162      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
163    register double u_i, u_min_i, u_max_i, u_j;    register double u_i, u_min_i, u_max_i, u_j;
164    register index_t iptr_ij;    register index_t iptr_ij;
   
   Paso_SystemMatrix_startCollect(L,u);  
165    n=Paso_SystemMatrix_getTotalNumRows(L);    n=Paso_SystemMatrix_getTotalNumRows(L);
166    pattern=L->pattern;    pattern=L->pattern;
   remote_u=Paso_SystemMatrix_finishCollect(L);  
167    #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)
168    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
169       u_i=u[i];       u_i=u[i];
# Line 163  void Paso_SolverFCT_setQs(const double* Line 174  void Paso_SolverFCT_setQs(const double*
174           j=pattern->mainPattern->index[iptr_ij];           j=pattern->mainPattern->index[iptr_ij];
175           u_j=u[j];           u_j=u[j];
176           u_min_i=MIN(u_min_i,u_j);           u_min_i=MIN(u_min_i,u_j);
177           u_max_i=MIN(u_max_i,u_j);           u_max_i=MAX(u_max_i,u_j);
178       }       }
179       #pragma ivdep       #pragma ivdep
180       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) {
181            j=pattern->couplePattern->index[iptr_ij];            j=pattern->col_couplePattern->index[iptr_ij];
182            u_j=remote_u[j];            u_j=remote_u[j];
183            u_min_i=MIN(u_min_i,u_j);            u_min_i=MIN(u_min_i,u_j);
184            u_max_i=MIN(u_max_i,u_j);            u_max_i=MAX(u_max_i,u_j);
185        }        }
186        QN[i]=u_min_i-u_i;        QN[i]=u_min_i-u_i;
187        QP[i]=u_max_i-u_i;        QP[i]=u_max_i-u_i;
# Line 179  void Paso_SolverFCT_setQs(const double* Line 190  void Paso_SolverFCT_setQs(const double*
190    
191  /*  /*
192   *   *
193   * 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])
194   *   *        
195   * m=fc->mass matrix   * m=fc->mass matrix
196   * 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)
197   */   */
198  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)
199  {  {
200    dim_t n, j, i;    dim_t n, j, i;
   double *remote_u;  
201    index_t iptr_ij;    index_t iptr_ij;
202    const double b2=-b;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
203    register double m_ij, ml_ij, k_ij, u_i;    const double *u_last= Paso_Coupler_borrowLocalData(fc->u_coupler);
204    Paso_SystemMatrixPattern *pattern;    const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
205    Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);    const double *remote_u_last=Paso_Coupler_borrowRemoteData(fc->u_coupler);
206      const double f1=- dt * (1.-fc->theta);
207      const double f2=  dt * fc->theta;
208      register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;
209      const Paso_SystemMatrixPattern *pattern=fc->iteration_matrix->pattern;
210    n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
211    pattern=fc->iteration_matrix->pattern;    if ( (ABS(f1) >0 ) ) {
212    remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);       if ( (ABS(f2) >0 ) ) {
213            #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)
214    if ( (ABS(a) >0 ) ) {          for (i = 0; i < n; ++i) {
215        if ( (ABS(b) >0 ) ) {             u_i=u[i];
216           #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];
217           for (i = 0; i < n; ++i) {             #pragma ivdep
218              u_i=u[i];             for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
219              #pragma ivdep                j=pattern->mainPattern->index[iptr_ij];
220              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];
221                  j=pattern->mainPattern->index[iptr_ij];                d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
222                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];                                                  fc->iteration_matrix->mainBlock->val[iptr_ij]);
223                  k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];                d_u=u[j]-u_i;
224                  ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];                d_u_last=u_last[j]-u_last_i;
225                  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;
226              }             }
227              #pragma ivdep             #pragma ivdep
228              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) {
229                  j=pattern->couplePattern->index[iptr_ij];                j=pattern->col_couplePattern->index[iptr_ij];
230                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];                m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
231                  k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];                d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
232                  ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];                                               fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
233                  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;
234              }                d_u_last=remote_u_last[j]-u_last_i;
235           }                flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
236        } else {             }
237           #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij)          }
238           for (i = 0; i < n; ++i) {       } else {
239              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)
240              #pragma ivdep          for (i = 0; i < n; ++i) {
241              for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {             u_i=u[i];
242                  j=pattern->mainPattern->index[iptr_ij];             u_last_i=u_last[i];
243                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];             #pragma ivdep
244                  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) {
245              }                j=pattern->mainPattern->index[iptr_ij];
246              #pragma ivdep                m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
247              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]+
248                  j=pattern->couplePattern->index[iptr_ij];                                                  fc->iteration_matrix->mainBlock->val[iptr_ij]);
249                  m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];                d_u=u[j]-u_i;
250                  flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*a*m_ij;                d_u_last=u_last[j]-u_last_i;
251              }                flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
252           }             }
253               #pragma ivdep
254               for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
255        }                j=pattern->col_couplePattern->index[iptr_ij];
256                  m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
257                  d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
258                                                 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
259                  d_u=remote_u[j]-u_i;
260                  d_u_last=remote_u_last[j]-u_last_i;
261                  flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
262               }
263            }
264         }
265    } else {    } else {
266        if ( (ABS(b) >0 ) ) {       if ( (ABS(f2) >0 ) ) {
267           #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)
268           for (i = 0; i < n; ++i) {          for (i = 0; i < n; ++i) {
269              u_i=u[i];             u_i=u[i];
270              #pragma ivdep             u_last_i=u_last[i];
271              for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {             #pragma ivdep
272                  j=pattern->mainPattern->index[iptr_ij];             for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
273                  k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];                j=pattern->mainPattern->index[iptr_ij];
274                  ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];                m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
275                  flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*b2*(ml_ij+k_ij);                d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
276              }                                                  fc->iteration_matrix->mainBlock->val[iptr_ij]);
277              #pragma ivdep                d_u=u[j]-u_i;
278              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;
279                  j=pattern->couplePattern->index[iptr_ij];                flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
280                  k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];             }
281                  ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];             #pragma ivdep
282                  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) {
283              }                j=pattern->col_couplePattern->index[iptr_ij];
284           }                m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
285                  d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
286        }                                               fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
287                  d_u=remote_u[j]-u_i;
288                  d_u_last=remote_u_last[j]-u_last_i;
289                  flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
290               }
291            }
292         } else {
293            #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)
294            for (i = 0; i < n; ++i) {
295               u_i=u[i];
296               u_last_i=u_last[i];
297               #pragma ivdep
298               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
299                  j=pattern->mainPattern->index[iptr_ij];
300                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
301                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
302                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
303                  d_u=u[j]-u_i;
304                  d_u_last=u_last[j]-u_last_i;
305                  flux_matrix->mainBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
306               }
307               #pragma ivdep
308               for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
309                  j=pattern->col_couplePattern->index[iptr_ij];
310                  m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
311                  d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
312                                                 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
313                  d_u=remote_u[j]-u_i;
314                  d_u_last=remote_u_last[j]-u_last_i;
315                  flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
316               }
317            }
318          }
319    }    }
320  }  }
321  /*  /*
322   *  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
323   *     *  
324   */   */
325  void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u)  void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const Paso_Coupler* u_coupler)
326  {  {
327    dim_t n, i, j;    dim_t n, i, j;
328    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
329    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
330      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
331    register double u_i, f_ij;    register double u_i, f_ij;
332    register index_t iptr_ij;    register index_t iptr_ij;
333    
334    n=Paso_SystemMatrix_getTotalNumRows(f);    n=Paso_SystemMatrix_getTotalNumRows(f);
335    pattern=f->pattern;    pattern=f->pattern;
   remote_u=Paso_SystemMatrix_finishCollect(f);  
336    #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)
337    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
338        u_i=u[i];        u_i=u[i];
# Line 288  void Paso_FCTransportProblem_applyPreAnt Line 343  void Paso_FCTransportProblem_applyPreAnt
343            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;
344        }        }
345        #pragma ivdep        #pragma ivdep
346        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) {
347            j=pattern->couplePattern->index[iptr_ij];            j=pattern->col_couplePattern->index[iptr_ij];
348            f_ij=f->coupleBlock->val[iptr_ij];            f_ij=f->col_coupleBlock->val[iptr_ij];
349            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;
350        }        }
351    }    }
352  }  }
353    
354    
355  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,
356                                     const double* QN,const double* QP,double* RN,double* RP)                                     const Paso_Coupler* QN_coupler,const Paso_Coupler* QP_coupler,double* RN,double* RP)
357  {  {
358    dim_t n, i, j;    dim_t n, i, j;
359    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
360      const double *QN=Paso_Coupler_borrowLocalData(QN_coupler);
361      const double *QP=Paso_Coupler_borrowLocalData(QP_coupler);
362    register double f_ij, PP_i, PN_i;    register double f_ij, PP_i, PN_i;
363    register index_t iptr_ij;    register index_t iptr_ij;
364    
# Line 316  void Paso_FCTransportProblem_setRs(const Line 373  void Paso_FCTransportProblem_setRs(const
373            j=pattern->mainPattern->index[iptr_ij];            j=pattern->mainPattern->index[iptr_ij];
374            if (i != j ) {            if (i != j ) {
375               f_ij=f->mainBlock->val[iptr_ij];               f_ij=f->mainBlock->val[iptr_ij];
376               if (f_ij <=0  ) {               if (f_ij <=0) {
377                  PN_i+=f_ij;                  PN_i+=f_ij;
378               } else {               } else {
379                  PP_i+=f_ij;                  PP_i+=f_ij;
# Line 324  void Paso_FCTransportProblem_setRs(const Line 381  void Paso_FCTransportProblem_setRs(const
381            }            }
382        }        }
383        #pragma ivdep        #pragma ivdep
384        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) {
385            j=pattern->couplePattern->index[iptr_ij];            j=pattern->col_couplePattern->index[iptr_ij];
386            f_ij=f->coupleBlock->val[iptr_ij];            f_ij=f->col_coupleBlock->val[iptr_ij];
387            if (f_ij <=0  ) {            if (f_ij <=0  ) {
388                  PN_i+=f_ij;                  PN_i+=f_ij;
389            } else {            } else {
# Line 347  void Paso_FCTransportProblem_setRs(const Line 404  void Paso_FCTransportProblem_setRs(const
404    
405  }  }
406    
407  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)
408  {  {
409    dim_t n, i, j;    dim_t n, i, j;
410    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
   double *remote_RN, *remote_RP;  
411    register double RN_i, RP_i, f_i, f_ij;    register double RN_i, RP_i, f_i, f_ij;
412    register index_t iptr_ij;    register index_t iptr_ij;
413      const double *RN=Paso_Coupler_borrowLocalData(RN_coupler);
414      const double *remote_RN=Paso_Coupler_borrowRemoteData(RN_coupler);
415      const double *RP=Paso_Coupler_borrowLocalData(RP_coupler);
416      const double *remote_RP=Paso_Coupler_borrowRemoteData(RP_coupler);
417    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
418    
   Paso_SystemMatrix_startCollect(flux_matrix,RN);  
419    pattern=flux_matrix->pattern;    pattern=flux_matrix->pattern;
   remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix);  
420    #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)
421    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
422       RN_i=RN[i];       RN_i=RN[i];
# Line 375  void Paso_FCTransportProblem_addCorrecte Line 433  void Paso_FCTransportProblem_addCorrecte
433           }           }
434       }       }
435       #pragma ivdep       #pragma ivdep
436       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) {
437            j=pattern->couplePattern->index[iptr_ij];            j=pattern->col_couplePattern->index[iptr_ij];
438            f_ij=flux_matrix->coupleBlock->val[iptr_ij];            f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];
439            if (f_ij >=0) {            if (f_ij >=0) {
440                f_i+=f_ij*MIN(RP_i,remote_RN[j]);                f_i+=f_ij*MIN(RP_i,remote_RN[j]);
441            }            }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) {  
442                f_i+=f_ij*MIN(RN_i,remote_RP[j]);                f_i+=f_ij*MIN(RN_i,remote_RP[j]);
443            }            }
444       }        }
445       f[i]+=f_i;        f[i]+=f_i;
446    }    }
447  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26