/[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 1371 by gross, Thu Jan 3 06:11:21 2008 UTC revision 1562 by gross, Wed May 21 13:04:40 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;    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              #pragma ivdep
67              for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
68                  j=pattern->mainPattern->index[iptr_ij];
69                  if (j!=i) {
70                     /* find entry a[j,i] */
71                     for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
72                        if (pattern->mainPattern->index[iptr_ji]==i) {
73                            rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];
74                            rtmp2=fc->transport_matrix->mainBlock->val[iptr_ji];
75                            d_ij=-MIN3(0.,rtmp1,rtmp2);
76                            fc->iteration_matrix->mainBlock->val[iptr_ij]=-(rtmp1+d_ij);
77                            sum-=d_ij;
78                            break;
79                        }
80                     }
81                 }
82             }
83              #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 */
100             fc->main_diagonal_low_order_transport_matrix[i]=sum;
101          }
102      }
103  }  }
104    /*
105  void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {   * out_i=m_i u_i + alpha \sum_{j <> i} l_{ij} (u_j-u_i) + beta q_i
106     *
107    double *remote_u=NULL;   */
108      void Paso_SolverFCT_setMuPaLuPbQ(double* out,
109    if (fc==NULL) return;                                   const double* M,
110    Paso_SystemMatrix_startCollect(fc->transport_matrix,u);                                   const Paso_Coupler* u_coupler,
111    /* process main block */                                   const double a,
112    Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);                                   const Paso_SystemMatrix *L,
113    /* finish exchange */                                   const double b,
114    remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);                                   const double* Q)
115    /* process couple block */  {
116    Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);    dim_t n, i, j;
117      Paso_SystemMatrixPattern *pattern;
118    Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);    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;
121      register index_t iptr_ij;
122    
123      n=Paso_SystemMatrix_getTotalNumRows(L);
124    
125      #pragma omp parallel for private(i) schedule(static)
126      for (i = 0; i < n; ++i) {
127          out[i]=M[i]*u[i]+b*Q[i];
128      }
129      if (ABS(a)>0) {
130          pattern=L->pattern;
131          #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
132          for (i = 0; i < n; ++i) {
133              sum=0;
134              u_i=u[i];
135              #pragma ivdep
136          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
137                   j=pattern->mainPattern->index[iptr_ij];
138                   l_ij=L->mainBlock->val[iptr_ij];
139                   sum+=l_ij*(u[j]-u_i);
140              }
141              #pragma ivdep
142          for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
143                   j=pattern->col_couplePattern->index[iptr_ij];
144                   l_ij=L->col_coupleBlock->val[iptr_ij];
145                   sum+=l_ij*(remote_u[j]-u_i);
146              }
147              out[i]+=a*sum;
148          }
149      }
150  }  }
151  /**************************************************************/  /*
152     * calculate QP[i] max_{i in L->pattern[i]} (u[j]-u[i] )
153  /* adds antidiffusion to fa     *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
154     *
    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  
   
155  */  */
156    void Paso_SolverFCT_setQs(const Paso_Coupler* u_coupler,double* QN, double* QP, const Paso_SystemMatrix *L)
157    {
158      dim_t n, i, j;
159      Paso_SystemMatrixPattern *pattern;
160      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;
163      register index_t iptr_ij;
164      n=Paso_SystemMatrix_getTotalNumRows(L);
165      pattern=L->pattern;
166      #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) {
168         u_i=u[i];
169         u_min_i=u_i;
170         u_max_i=u_i;
171         #pragma ivdep
172         for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
173             j=pattern->mainPattern->index[iptr_ij];
174             u_j=u[j];
175             u_min_i=MIN(u_min_i,u_j);
176             u_max_i=MAX(u_max_i,u_j);
177         }
178         #pragma ivdep
179         for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
180              j=pattern->col_couplePattern->index[iptr_ij];
181              u_j=remote_u[j];
182              u_min_i=MIN(u_min_i,u_j);
183              u_max_i=MAX(u_max_i,u_j);
184          }
185          QN[i]=u_min_i-u_i;
186          QP[i]=u_max_i-u_i;
187      }
188    }
189    
190  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {  /*
191     *
192    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;   *  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    index_t color, iptr_ij,j,iptr_ji, i;   *        
194    dim_t n;   * m=fc->mass matrix
195     * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
196     */
197    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;
200      index_t iptr_ij;
201      const double *u=Paso_Coupler_borrowLocalData(u_coupler);
202      const double *u_last= Paso_Coupler_borrowLocalData(fc->u_coupler);
203      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
204      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);
210      if ( (ABS(f1) >0 ) ) {
211         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            for (i = 0; i < n; ++i) {
214               u_i=u[i];
215               u_last_i=u_last[i];
216               #pragma ivdep
217               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
218                  j=pattern->mainPattern->index[iptr_ij];
219                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
220                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
221                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
222                  d_u=u[j]-u_i;
223                  d_u_last=u_last[j]-u_last_i;
224                  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
227               for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
228                  j=pattern->col_couplePattern->index[iptr_ij];
229                  m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
230                  d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
231                                                 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
232                  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               }
236            }
237         } else {
238            #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            for (i = 0; i < n; ++i) {
240               u_i=u[i];
241               u_last_i=u_last[i];
242               #pragma ivdep
243               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
244                  j=pattern->mainPattern->index[iptr_ij];
245                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
246                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
247                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
248                  d_u=u[j]-u_i;
249                  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 {
265         if ( (ABS(f2) >0 ) ) {
266            #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) {
268               u_i=u[i];
269               u_last_i=u_last[i];
270               #pragma ivdep
271               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
272                  j=pattern->mainPattern->index[iptr_ij];
273                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
274                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
275                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
276                  d_u=u[j]-u_i;
277                  d_u_last=u_last[j]-u_last_i;
278                  flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
279               }
280               #pragma ivdep
281               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
322     *  
323     */
324    void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const Paso_Coupler* u_coupler)
325    {
326      dim_t n, i, j;
327      Paso_SystemMatrixPattern *pattern;
328      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;
331      register index_t iptr_ij;
332    
333      n=Paso_SystemMatrix_getTotalNumRows(f);
334      pattern=f->pattern;
335      #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
336      for (i = 0; i < n; ++i) {
337          u_i=u[i];
338          #pragma ivdep
339          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
340              j=pattern->mainPattern->index[iptr_ij];
341              f_ij=f->mainBlock->val[iptr_ij];
342              if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
343          }
344          #pragma ivdep
345          for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
346              j=pattern->col_couplePattern->index[iptr_ij];
347              f_ij=f->col_coupleBlock->val[iptr_ij];
348              if (f_ij * (u_i-remote_u[j]) <= 0) f->col_coupleBlock->val[iptr_ij]=0;
349          }
350      }
351    }
352    
   if (fc==NULL) return;  
   n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);  
   /* exchange */  
353    
354    void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
355                                       const Paso_Coupler* QN_coupler,const Paso_Coupler* QP_coupler,double* RN,double* RP)
356    {
357      dim_t n, i, j;
358      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;
362      register index_t iptr_ij;
363    
364      n=Paso_SystemMatrix_getTotalNumRows(f);
365      pattern=f->pattern;
366      #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
367      for (i = 0; i < n; ++i) {
368          PP_i=0.;
369          PN_i=0.;
370          #pragma ivdep
371          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
372              j=pattern->mainPattern->index[iptr_ij];
373              if (i != j ) {
374                 f_ij=f->mainBlock->val[iptr_ij];
375                 if (f_ij <=0) {
376                    PN_i+=f_ij;
377                 } else {
378                    PP_i+=f_ij;
379                 }
380              }
381          }
382          #pragma ivdep
383          for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
384              j=pattern->col_couplePattern->index[iptr_ij];
385              f_ij=f->col_coupleBlock->val[iptr_ij];
386              if (f_ij <=0  ) {
387                    PN_i+=f_ij;
388              } else {
389                    PP_i+=f_ij;
390              }
391          }
392          if (PN_i<0) {
393             RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i);
394          } else {
395             RN[i]=1.;
396          }
397          if (PP_i>0) {
398             RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i);
399          } else {
400             RP[i]=1.;
401          }
402      }
403    
404    #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, d_ij)  
            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.;  
                   r_p=0.;  
                   r_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;  
 printf("%d %d : %e %e :: %e \n",i,j,u_i,u[j],a_ij);  
                       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.;  
 printf("%d: %e %e %e : %e %e %e : %e\n",i,Q_p,P_p,r_p,Q_n,P_n,r_n,u_i);  
                   /* 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.) && ( (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*d_ij,a_ji+d_ij)*(u_i-u_j);  
   
                                     fa[i]+=f_ij;  
                                     fa[j]-=f_ij;  
 printf("%d %d => %e %e : %e %e : %e %e : fa[%d]=%e fa[%d]=%e\n",i,j,d_ij,(u_i-u_j), a_ij, a_ji, r_ij,f_ij,i,fa[i],j,fa[j]);  
   
                                     
                                 }  
                                 break;  
                             }  
                         }  
                      }  
                   }  
                   /* anti diffusive flux from couple block */  
405    
406                    /* TODO */  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;
409             #pragma omp barrier    Paso_SystemMatrixPattern *pattern;
410         }    register double RN_i, RP_i, f_i, f_ij;
411      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);
417    
418      pattern=flux_matrix->pattern;
419      #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) {
421         RN_i=RN[i];
422         RP_i=RP[i];
423         f_i=0;
424         #pragma ivdep
425         for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
426             j=pattern->mainPattern->index[iptr_ij];
427             f_ij=flux_matrix->mainBlock->val[iptr_ij];
428             if (f_ij >=0) {
429                  f_i+=f_ij*MIN(RP_i,RN[j]);
430             } else {
431                  f_i+=f_ij*MIN(RN_i,RP[j]);
432             }
433         }
434         #pragma ivdep
435         for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
436              j=pattern->col_couplePattern->index[iptr_ij];
437              f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];
438              if (f_ij >=0) {
439                  f_i+=f_ij*MIN(RP_i,remote_RN[j]);
440              }else {
441                  f_i+=f_ij*MIN(RN_i,remote_RP[j]);
442              }
443          }
444          f[i]+=f_i;
445    }    }
446  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26