/[escript]/trunk/paso/src/SolverFCT_FluxControl.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1370 by gross, Wed Jan 2 09:21:43 2008 UTC revision 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 26  Line 27 
27  #include "SolverFCT.h"  #include "SolverFCT.h"
28  #include "PasoUtil.h"  #include "PasoUtil.h"
29    
 #define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.)  
 #define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b)))  
 #define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b))))  
   
 #define FLUX_L(a,b) SUPERBEE(a,b)  /* alter for other flux limiter */  
   
 #define FLUX_LIMITER(a) FLUX_L(a,1)  
30    
31  /**************************************************************/  /**************************************************************/
32    
33  /* adds A plus stabelising diffusion into the matrix B        */  /* create the low order transport matrix and stores it negative values
34     * into the iteration_matrix accept the main diagonal which is stored seperately
35  /* d_ij=alpha*max(0,-a[i,j],-a[j,i])  */   * if fc->iteration_matrix==NULL, fc->iteration_matrix is allocated
36  /* b[i,j]+=alpha*(a[i,j]+d_ij)  */   *
37  /* b[j,i]+=alpha*(a[j,i]+d_ij)  */   * a=transport_matrix
38  /* b[i,i]-=alpha*d_ij  */   * b= low_order_transport_matrix = - iteration_matrix
39  /* b[j,j]-=alpha*d_ij  */   * c=main diagonal low_order_transport_matrix
40     * initialize c[i] mit a[i,i]
41  void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {   *    d_ij=max(0,-a[i,j],-a[j,i])  
42    dim_t n,i;   *    b[i,j]=-(a[i,j]+d_ij)        
43    index_t color, iptr_ij,j,iptr_ji;   *    c[i]-=d_ij                  
44    register double d_ij;   */
45    
46    void Paso_FCTransportProblem_setLowOrderOperator(Paso_FCTransportProblem * fc) {
47      dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i;
48      index_t iptr_ij,j,iptr_ji;
49      Paso_SystemMatrixPattern *pattern;
50      register double d_ij, sum, rtmp1, rtmp2;
51    
52    if (fc==NULL) return;    if (fc==NULL) return;
53    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);    if (fc->iteration_matrix==NULL) {
54          fc->iteration_matrix=Paso_SystemMatrix_alloc(fc->transport_matrix->type,
55    #pragma omp parallel private(color)                                                     fc->transport_matrix->pattern,
56    {                                                     fc->transport_matrix->row_block_size,
57         /* 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  
        }  
58    }    }
59    
60      if (Paso_noError()) {
61          pattern=fc->iteration_matrix->pattern;
62          n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
63          #pragma omp parallel for private(i,sum,iptr_ij,j,iptr_ji,rtmp1, rtmp2,d_ij)  schedule(static)
64          for (i = 0; i < n; ++i) {
65              sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
66              /* 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) {
69                  j=pattern->mainPattern->index[iptr_ij];
70                  if (j!=i) {
71                     /* find entry a[j,i] */
72                     for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
73                        if (pattern->mainPattern->index[iptr_ji]==i) {
74                            rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];
75                            rtmp2=fc->transport_matrix->mainBlock->val[iptr_ji];
76                            d_ij=-MIN3(0.,rtmp1,rtmp2);
77                            fc->iteration_matrix->mainBlock->val[iptr_ij]=-(rtmp1+d_ij);
78                            sum-=d_ij;
79                            break;
80                        }
81                     }
82                 }
83             }
84              #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 */
101             fc->main_diagonal_low_order_transport_matrix[i]=sum;
102          }
103      }
104  }  }
105    /*
106  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
107     *
108    double *remote_u=NULL;   */
109      void Paso_SolverFCT_setMuPaLuPbQ(double* out,
110    if (fc==NULL) return;                                   const double* M,
111    Paso_SystemMatrix_startCollect(fc->transport_matrix,u);                                   const Paso_Coupler* u_coupler,
112    /* process main block */                                   const double a,
113    Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);                                   const Paso_SystemMatrix *L,
114    /* finish exchange */                                   const double b,
115    remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);                                   const double* Q)
116    /* process couple block */  {
117    Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);    dim_t n, i, j;
118      Paso_SystemMatrixPattern *pattern;
119    Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);    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;
122      register index_t iptr_ij;
123    
124      n=Paso_SystemMatrix_getTotalNumRows(L);
125    
126      #pragma omp parallel for private(i) schedule(static)
127      for (i = 0; i < n; ++i) {
128          out[i]=M[i]*u[i]+b*Q[i];
129      }
130      if (ABS(a)>0) {
131          pattern=L->pattern;
132          #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
133          for (i = 0; i < n; ++i) {
134              sum=0;
135              u_i=u[i];
136              #pragma ivdep
137          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
138                   j=pattern->mainPattern->index[iptr_ij];
139                   l_ij=L->mainBlock->val[iptr_ij];
140                   sum+=l_ij*(u[j]-u_i);
141              }
142              #pragma ivdep
143          for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
144                   j=pattern->col_couplePattern->index[iptr_ij];
145                   l_ij=L->col_coupleBlock->val[iptr_ij];
146                   sum+=l_ij*(remote_u[j]-u_i);
147              }
148              out[i]+=a*sum;
149          }
150      }
151  }  }
152  /**************************************************************/  /*
153     * calculate QP[i] max_{i in L->pattern[i]} (u[j]-u[i] )
154  /* adds antidiffusion to fa     *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
155     *
    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  
   
156  */  */
157    void Paso_SolverFCT_setQs(const Paso_Coupler* u_coupler,double* QN, double* QP, const Paso_SystemMatrix *L)
158    {
159      dim_t n, i, j;
160      Paso_SystemMatrixPattern *pattern;
161      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;
164      register index_t iptr_ij;
165      n=Paso_SystemMatrix_getTotalNumRows(L);
166      pattern=L->pattern;
167      #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) {
169         u_i=u[i];
170         u_min_i=u_i;
171         u_max_i=u_i;
172         #pragma ivdep
173         for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
174             j=pattern->mainPattern->index[iptr_ij];
175             u_j=u[j];
176             u_min_i=MIN(u_min_i,u_j);
177             u_max_i=MAX(u_max_i,u_j);
178         }
179         #pragma ivdep
180         for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
181              j=pattern->col_couplePattern->index[iptr_ij];
182              u_j=remote_u[j];
183              u_min_i=MIN(u_min_i,u_j);
184              u_max_i=MAX(u_max_i,u_j);
185          }
186          QN[i]=u_min_i-u_i;
187          QP[i]=u_max_i-u_i;
188      }
189    }
190    
191  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {  /*
192     *
193    register double u_i,P_p,P_n,Q_p,Q_n,r_p,r_n, a_ij, d, u_j, r_ij, f_ij, a_ji;   *  f_{ij} = (m_{ij} - dt (1-theta) d_{ij}) (u_last[j]-u_last[i]) - (m_{ij} + dt theta d_{ij}) (u[j]-u[i])
194    index_t color, iptr_ij,j,iptr_ji, i;   *        
195    dim_t n;   * m=fc->mass matrix
196     * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
197     */
198    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;
201      index_t iptr_ij;
202      const double *u=Paso_Coupler_borrowLocalData(u_coupler);
203      const double *u_last= Paso_Coupler_borrowLocalData(fc->u_coupler);
204      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
205      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);
211      if ( (ABS(f1) >0 ) ) {
212         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            for (i = 0; i < n; ++i) {
215               u_i=u[i];
216               u_last_i=u_last[i];
217               #pragma ivdep
218               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
219                  j=pattern->mainPattern->index[iptr_ij];
220                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
221                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
222                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
223                  d_u=u[j]-u_i;
224                  d_u_last=u_last[j]-u_last_i;
225                  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
228               for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
229                  j=pattern->col_couplePattern->index[iptr_ij];
230                  m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
231                  d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
232                                                 fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
233                  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               }
237            }
238         } else {
239            #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            for (i = 0; i < n; ++i) {
241               u_i=u[i];
242               u_last_i=u_last[i];
243               #pragma ivdep
244               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
245                  j=pattern->mainPattern->index[iptr_ij];
246                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
247                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
248                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
249                  d_u=u[j]-u_i;
250                  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 {
266         if ( (ABS(f2) >0 ) ) {
267            #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) {
269               u_i=u[i];
270               u_last_i=u_last[i];
271               #pragma ivdep
272               for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
273                  j=pattern->mainPattern->index[iptr_ij];
274                  m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
275                  d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
276                                                    fc->iteration_matrix->mainBlock->val[iptr_ij]);
277                  d_u=u[j]-u_i;
278                  d_u_last=u_last[j]-u_last_i;
279                  flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
280               }
281               #pragma ivdep
282               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
323     *  
324     */
325    void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const Paso_Coupler* u_coupler)
326    {
327      dim_t n, i, j;
328      Paso_SystemMatrixPattern *pattern;
329      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;
332      register index_t iptr_ij;
333    
334      n=Paso_SystemMatrix_getTotalNumRows(f);
335      pattern=f->pattern;
336      #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
337      for (i = 0; i < n; ++i) {
338          u_i=u[i];
339          #pragma ivdep
340          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
341              j=pattern->mainPattern->index[iptr_ij];
342              f_ij=f->mainBlock->val[iptr_ij];
343              if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
344          }
345          #pragma ivdep
346          for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
347              j=pattern->col_couplePattern->index[iptr_ij];
348              f_ij=f->col_coupleBlock->val[iptr_ij];
349              if (f_ij * (u_i-remote_u[j]) <= 0) f->col_coupleBlock->val[iptr_ij]=0;
350          }
351      }
352    }
353    
   if (fc==NULL) return;  
   n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);  
   /* exchange */  
354    
355    void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
356                                       const Paso_Coupler* QN_coupler,const Paso_Coupler* QP_coupler,double* RN,double* RP)
357    {
358      dim_t n, i, j;
359      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;
363      register index_t iptr_ij;
364    
365      n=Paso_SystemMatrix_getTotalNumRows(f);
366      pattern=f->pattern;
367      #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
368      for (i = 0; i < n; ++i) {
369          PP_i=0.;
370          PN_i=0.;
371          #pragma ivdep
372          for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
373              j=pattern->mainPattern->index[iptr_ij];
374              if (i != j ) {
375                 f_ij=f->mainBlock->val[iptr_ij];
376                 if (f_ij <=0) {
377                    PN_i+=f_ij;
378                 } else {
379                    PP_i+=f_ij;
380                 }
381              }
382          }
383          #pragma ivdep
384          for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
385              j=pattern->col_couplePattern->index[iptr_ij];
386              f_ij=f->col_coupleBlock->val[iptr_ij];
387              if (f_ij <=0  ) {
388                    PN_i+=f_ij;
389              } else {
390                    PP_i+=f_ij;
391              }
392          }
393          if (PN_i<0) {
394             RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i);
395          } else {
396             RN[i]=1.;
397          }
398          if (PP_i>0) {
399             RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i);
400          } else {
401             RP[i]=1.;
402          }
403      }
404    
405    #pragma omp parallel private(color)  }
   {  
        for (color=0;color<fc->num_colors;++color) {  
            #pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,r_p,r_n,iptr_ij,a_ij,d,j,iptr_ji, u_j, r_ij, f_ij, a_ji)  
            for (i = 0; i < n; ++i) {  
               if (fc->colorOf[i]==color) {  
                   u_i=u[i];  
                   /* gather the smoothness sensor */  
                   P_p=0.;  
                   P_n=0.;  
                   Q_p=0.;  
                   Q_n=0.;  
                   #pragma ivdep  
               for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {  
                       a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];  
                       j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];  
                       d=u[j]-u_i;  
                       if (a_ij<0.) {  
                          if (d<0.) {  
                             P_p+=a_ij*d;  
                          } else {  
                             P_n+=a_ij*d;  
                          }  
                       } else {  
                          if (d>0.) {  
                             Q_p+=a_ij*d;  
                          } else {  
                             Q_n+=a_ij*d;  
                          }  
                       }  
               }  
                   #pragma ivdep  
               for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {  
                       a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];  
                       j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];  
                       d=u_remote[j]-u_i;  
                       if (a_ij<0.) {  
                          if (d<0.) {  
                             P_p+=a_ij*d;  
                          } else {  
                             P_n+=a_ij*d;  
                          }  
                       } else {  
                          if (d>0.) {  
                             Q_p+=a_ij*d;  
                          } else {  
                             Q_n+=a_ij*d;  
                          }  
                       }  
               }  
                   /* set the smoothness indicators */  
                   r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;  
                   r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;  
                   /* anti diffusive flux from main block */  
                   for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {  
                      a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];  
                      j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];  
                      if (a_ij < 0 && i!=j) {  
                         /* find entry a[j,i] */  
                         for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {  
                             if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {  
                                 a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];  
                                 if  (a_ji > a_ij || (a_ji == a_ij && j<i) ) {  
                                     u_j=u[j];  
                                     r_ij = u_i>u_j ? r_p : r_n;  
                                     f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j);  
                                     fa[i]+=f_ij;  
                                     fa[j]-=f_ij;  
                                     break;  
                                 }  
                             }  
                         }  
                      }  
                   }  
                   /* anti diffusive flux from couple block */  
406    
407                    /* TODO */  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;
410             #pragma omp barrier    Paso_SystemMatrixPattern *pattern;
411         }    register double RN_i, RP_i, f_i, f_ij;
412      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);
418    
419      pattern=flux_matrix->pattern;
420      #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) {
422         RN_i=RN[i];
423         RP_i=RP[i];
424         f_i=0;
425         #pragma ivdep
426         for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
427             j=pattern->mainPattern->index[iptr_ij];
428             f_ij=flux_matrix->mainBlock->val[iptr_ij];
429             if (f_ij >=0) {
430                  f_i+=f_ij*MIN(RP_i,RN[j]);
431             } else {
432                  f_i+=f_ij*MIN(RN_i,RP[j]);
433             }
434         }
435         #pragma ivdep
436         for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
437              j=pattern->col_couplePattern->index[iptr_ij];
438              f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];
439              if (f_ij >=0) {
440                  f_i+=f_ij*MIN(RP_i,remote_RN[j]);
441              }else {
442                  f_i+=f_ij*MIN(RN_i,remote_RP[j]);
443              }
444          }
445          f[i]+=f_i;
446    }    }
447  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26