/[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 1401 by gross, Fri Jan 25 04:31:18 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 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)  
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, sum;   */
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                                                       fc->transport_matrix->pattern,
56                                                       fc->transport_matrix->row_block_size,
57                                                       fc->transport_matrix->col_block_size);
58      }
59    
60    #pragma omp parallel for private(i,iptr_ij,j,iptr_ji,d_ij,sum)  schedule(static)    if (Paso_noError()) {
61    for (i = 0; i < n; ++i) {        pattern=fc->iteration_matrix->pattern;
62       sum=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];        n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
63       for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {        #pragma omp parallel for private(i,sum,iptr_ij,j,iptr_ji,rtmp1, rtmp2,d_ij)  schedule(static)
64           j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];        for (i = 0; i < n; ++i) {
65           if (j!=i) {            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
66               /* find entry a[j,i] */            /* look at a[i,j] */
67               for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {            #pragma ivdep
68                  if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
69                      d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],                j=pattern->mainPattern->index[iptr_ij];
70                                            fc->flux_matrix->mainBlock->val[iptr_ji]);                if (j!=i) {
71                      fc->transport_matrix->mainBlock->val[iptr_ij]+=                   /* find entry a[j,i] */
72                                           alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;                   for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
73                      sum-=d_ij;                      if (pattern->mainPattern->index[iptr_ji]==i) {
74                      break;                          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       /* TODO process couple block */            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       /* update main diagonal */                /* find entry a[j,i] */
88       fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=sum;                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     * out_i=m_i u_i + alpha \sum_{j <> i} l_{ij} (u_j-u_i) + beta q_i
107     *
108     */
109    void Paso_SolverFCT_setMuPaLuPbQ(double* out,
110                                     const double* M,
111                                     const Paso_Coupler* u_coupler,
112                                     const double a,
113                                     const Paso_SystemMatrix *L,
114                                     const double b,
115                                     const double* Q)
116    {
117      dim_t n, i, j;
118      Paso_SystemMatrixPattern *pattern;
119      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  void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {    n=Paso_SystemMatrix_getTotalNumRows(L);
   /*  
    *   sets fa=transport_matrix*u+anti-diffuison_flux(u)  
    */  
   
   double *remote_u=NULL;  
     
   if (fc==NULL) return;  
   Paso_SystemMatrix_startCollect(fc->transport_matrix,u);  
   /* process main block */  
   Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);  
   /* finish exchange */  
   remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);  
   /* process couple block */  
   Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);  
125    
126    Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);    #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, d_ij, sum;   *  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    if (fc==NULL) return;  void Paso_FCTransportProblem_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const Paso_Coupler* u_coupler)
199    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);  {
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    
354    
355    #pragma omp parallel  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         * calculate the smootness sensors    dim_t n, i, j;
359        */    Paso_SystemMatrixPattern *pattern;
360        #pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,iptr_ij,a_ij,j,d)    const double *QN=Paso_Coupler_borrowLocalData(QN_coupler);
361        for (i = 0; i < n; ++i) {    const double *QP=Paso_Coupler_borrowLocalData(QP_coupler);
362            u_i=u[i];    register double f_ij, PP_i, PN_i;
363            P_p=0.;    register index_t iptr_ij;
364            P_n=0.;  
365            Q_p=0.;    n=Paso_SystemMatrix_getTotalNumRows(f);
366            Q_n=0.;    pattern=f->pattern;
367            #pragma ivdep    #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
368        for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {    for (i = 0; i < n; ++i) {
369                 a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];        PP_i=0.;
370                 j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];        PN_i=0.;
371                 d=u[j]-u_i;        #pragma ivdep
372                 if (a_ij<0.) {        for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
373                    if (d<0.) {            j=pattern->mainPattern->index[iptr_ij];
374                        P_p+=a_ij*d;            if (i != j ) {
375                    } else {               f_ij=f->mainBlock->val[iptr_ij];
376                        P_n+=a_ij*d;               if (f_ij <=0) {
377                    }                  PN_i+=f_ij;
378                 } else {               } else {
379                    if (d>0.) {                  PP_i+=f_ij;
380                      Q_p+=a_ij*d;               }
                   } else {  
                     Q_n+=a_ij*d;  
                   }  
                }  
       }  
           #pragma ivdep  
       for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {  
                a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];  
                j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];  
                d=u_remote[j]-u_i;  
                if (a_ij<0.) {  
                    if (d<0.) {  
                      P_p+=a_ij*d;  
                    } else {  
                      P_n+=a_ij*d;  
                    }  
                } else {  
                   if (d>0.) {  
                      Q_p+=a_ij*d;  
                   } else {  
                      Q_n+=a_ij*d;  
                   }  
                }  
       }  
           /* set the smoothness indicators */  
           fc->r_p[i] = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;  
           fc->r_n[i] = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.;  
   
       } /* end of row loop for smootheness indicator */  
   
       /*  
        * calculate antidiffusion  
       */  
       #pragma omp for schedule(static) private(i, u_i, sum, iptr_ij, a_ij, j, iptr_ji, a_ji,d_ij, u_j, r_ij, f_ij)  
       for (i = 0; i < n; ++i) {  
           u_i=u[i];  
           sum=0;  
           /* anti diffusive flux from main block */  
           for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {  
               a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];  
               j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];  
               if ( i!=j ) {  
                    /* find entry a[j,i] */  
                    for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {  
                       if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {  
                            a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];  
                            d_ij=-MIN3(0,a_ji,a_ij);  
                            if (d_ij > 0.) {  
                                 u_j=u[j];  
                                 if (a_ji >= a_ij) {  
                                     r_ij = u_i>u_j ? fc->r_p[i] : fc->r_n[i];  
                                     f_ij =MIN(r_ij*d_ij,a_ji+d_ij);  
                                 } else {  
                                     r_ij = u_j>u_i ? fc->r_p[j] : fc->r_n[j];  
                                     f_ij =MIN(r_ij*d_ij,a_ij+d_ij);  
                                 }  
                                 sum+=f_ij*(u_i-u_j);  
                            }  
                            break;  
                       }  
                    }  
               }  
381            }            }
382            /* anti diffusive flux from couple block */        }
383            /* TODO */        #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    }
406    
407    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      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            fa[i]+=sum;    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    } /* end of parallel block */        f[i]+=f_i;
446      }
447  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26