/[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 1552 by gross, Thu May 8 08:52:41 2008 UTC revision 1562 by gross, Wed May 21 13:04:40 2008 UTC
# Line 63  void Paso_FCTransportProblem_setLowOrder Line 63  void Paso_FCTransportProblem_setLowOrder
63        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
64            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
65            /* look at a[i,j] */            /* look at a[i,j] */
66              #pragma ivdep
67            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
68                j=pattern->mainPattern->index[iptr_ij];                j=pattern->mainPattern->index[iptr_ij];
69                if (j!=i) {                if (j!=i) {
# Line 79  void Paso_FCTransportProblem_setLowOrder Line 80  void Paso_FCTransportProblem_setLowOrder
80                   }                   }
81               }               }
82           }           }
83           /* TODO process couple block */            #pragma ivdep
84                      for (iptr_ij=pattern->col_couplePattern->ptr[i];iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
85                  j=pattern->col_couplePattern->index[iptr_ij];
86                  /* find entry a[j,i] */
87                  for (iptr_ji=pattern->row_couplePattern->ptr[j]; iptr_ji<pattern->row_couplePattern->ptr[j+1]; ++iptr_ji) {
88                        if (pattern->row_couplePattern->index[iptr_ji]==i) {
89                            rtmp1=fc->transport_matrix->col_coupleBlock->val[iptr_ij];
90                            rtmp2=fc->transport_matrix->row_coupleBlock->val[iptr_ji];
91                            d_ij=-MIN3(0.,rtmp1,rtmp2);
92                            fc->iteration_matrix->col_coupleBlock->val[iptr_ij]=-(rtmp1+d_ij);
93                            fc->iteration_matrix->row_coupleBlock->val[iptr_ji]=-(rtmp2+d_ij);
94                            sum-=d_ij;
95                            break;
96                        }
97                  }
98             }
99           /* set main diagonal entry */           /* set main diagonal entry */
100           fc->main_diagonal_low_order_transport_matrix[i]=sum;           fc->main_diagonal_low_order_transport_matrix[i]=sum;
101        }        }
# Line 92  void Paso_FCTransportProblem_setLowOrder Line 107  void Paso_FCTransportProblem_setLowOrder
107   */   */
108  void Paso_SolverFCT_setMuPaLuPbQ(double* out,  void Paso_SolverFCT_setMuPaLuPbQ(double* out,
109                                   const double* M,                                   const double* M,
110                                   const double* u,                                   const Paso_Coupler* u_coupler,
111                                   const double a,                                   const double a,
112                                   Paso_SystemMatrix *L,                                   const Paso_SystemMatrix *L,
113                                   const double b,                                   const double b,
114                                   const double* Q)                                   const double* Q)
115  {  {
116    dim_t n, i, j;    dim_t n, i, j;
117    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
118    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
119      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
120    register double sum, u_i, l_ij;    register double sum, u_i, l_ij;
121    register index_t iptr_ij;    register index_t iptr_ij;
122    
123    n=Paso_SystemMatrix_getTotalNumRows(L);    n=Paso_SystemMatrix_getTotalNumRows(L);
124    
   if (ABS(a)>0) {  
       Paso_SystemMatrix_startCollect(L,u);  
   }  
125    #pragma omp parallel for private(i) schedule(static)    #pragma omp parallel for private(i) schedule(static)
126    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
127        out[i]=M[i]*u[i]+b*Q[i];        out[i]=M[i]*u[i]+b*Q[i];
128    }    }
129    if (ABS(a)>0) {    if (ABS(a)>0) {
130        pattern=L->pattern;        pattern=L->pattern;
       remote_u=Paso_SystemMatrix_finishCollect(L);  
131        #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)        #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
132        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
133            sum=0;            sum=0;
# Line 141  void Paso_SolverFCT_setMuPaLuPbQ(double* Line 153  void Paso_SolverFCT_setMuPaLuPbQ(double*
153   *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )   *           QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
154   *   *
155  */  */
156  void Paso_SolverFCT_setQs(const double* u,double* QN, double* QP, Paso_SystemMatrix *L)  void Paso_SolverFCT_setQs(const Paso_Coupler* u_coupler,double* QN, double* QP, const Paso_SystemMatrix *L)
157  {  {
158    dim_t n, i, j;    dim_t n, i, j;
159    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
160    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
161      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
162    register double u_i, u_min_i, u_max_i, u_j;    register double u_i, u_min_i, u_max_i, u_j;
163    register index_t iptr_ij;    register index_t iptr_ij;
   
   Paso_SystemMatrix_startCollect(L,u);  
164    n=Paso_SystemMatrix_getTotalNumRows(L);    n=Paso_SystemMatrix_getTotalNumRows(L);
165    pattern=L->pattern;    pattern=L->pattern;
   remote_u=Paso_SystemMatrix_finishCollect(L);  
166    #pragma omp parallel for schedule(static) private(i, u_i, u_min_i, u_max_i, j, u_j, iptr_ij)    #pragma omp parallel for schedule(static) private(i, u_i, u_min_i, u_max_i, j, u_j, iptr_ij)
167    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
168       u_i=u[i];       u_i=u[i];
# Line 184  void Paso_SolverFCT_setQs(const double* Line 194  void Paso_SolverFCT_setQs(const double*
194   * m=fc->mass matrix   * m=fc->mass matrix
195   * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)   * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
196   */   */
197  void Paso_FCTransportProblem_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const double* u, const double* u_last)  void Paso_FCTransportProblem_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const Paso_Coupler* u_coupler)
198  {  {
199    dim_t n, j, i;    dim_t n, j, i;
   double *remote_u=NULL, *remote_u_last=NULL;  
200    index_t iptr_ij;    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);    const double f1=- dt * (1.-fc->theta);
206    const double f2=  dt * fc->theta;    const double f2=  dt * fc->theta;
207    register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;    register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;
208    Paso_SystemMatrixPattern *pattern;    const Paso_SystemMatrixPattern *pattern=fc->iteration_matrix->pattern;
   Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);  
   Paso_SystemMatrix_startCollect(fc->iteration_matrix,u_last);  
209    n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
   pattern=fc->iteration_matrix->pattern;  
   remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);  
   remote_u_last=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);  
210    if ( (ABS(f1) >0 ) ) {    if ( (ABS(f1) >0 ) ) {
211       if ( (ABS(f2) >0 ) ) {       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)          #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)
# Line 310  void Paso_FCTransportProblem_setAntiDiff Line 318  void Paso_FCTransportProblem_setAntiDiff
318    }    }
319  }  }
320  /*  /*
  *  
  * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])  
  *  
  * m=fc->mass matrix  
  * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)  
  */  
 void Paso_FCTransportProblem_updateAntiDiffusionFlux(const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix,const double a, const double b, const double* u)  
 {  
   dim_t n, j, i;  
   double *remote_u;  
   index_t iptr_ij;  
   const double b2=-b;  
   register double m_ij, ml_ij, k_ij, u_i;  
   Paso_SystemMatrixPattern *pattern;  
   Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);  
   n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);  
   pattern=fc->iteration_matrix->pattern;  
   remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);  
   if ( (ABS(a) >0 ) ) {  
       if ( (ABS(b) >0 ) ) {  
          #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij)  
          for (i = 0; i < n; ++i) {  
             u_i=u[i];  
             #pragma ivdep  
             for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {  
                 j=pattern->mainPattern->index[iptr_ij];  
                 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];  
                 k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];  
                 ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];  
                 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));  
             }  
             #pragma ivdep  
             for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {  
                 j=pattern->col_couplePattern->index[iptr_ij];  
                 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];  
                 k_ij=fc->transport_matrix->col_coupleBlock->val[iptr_ij];  
                 ml_ij=fc->iteration_matrix->col_coupleBlock->val[iptr_ij];  
                 flux_matrix->col_coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));  
             }  
          }  
       } else {  
          #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij)  
          for (i = 0; i < n; ++i) {  
             u_i=u[i];  
             #pragma ivdep  
             for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {  
                 j=pattern->mainPattern->index[iptr_ij];  
                 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];  
                 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*a*m_ij;  
             }  
             #pragma ivdep  
             for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {  
                 j=pattern->col_couplePattern->index[iptr_ij];  
                 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];  
                 flux_matrix->col_coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*a*m_ij;  
             }  
          }  
   
   
       }  
   } else {  
       if ( (ABS(b) >0 ) ) {  
          #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,k_ij, ml_ij)  
          for (i = 0; i < n; ++i) {  
             u_i=u[i];  
             #pragma ivdep  
             for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {  
                 j=pattern->mainPattern->index[iptr_ij];  
                 k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];  
                 ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];  
                 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*b2*(ml_ij+k_ij);  
             }  
             #pragma ivdep  
             for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {  
                 j=pattern->col_couplePattern->index[iptr_ij];  
                 k_ij=fc->transport_matrix->col_coupleBlock->val[iptr_ij];  
                 ml_ij=fc->iteration_matrix->col_coupleBlock->val[iptr_ij];  
                 flux_matrix->col_coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*b2*(ml_ij+k_ij);  
             }  
          }  
   
       }  
   }  
 }  
 /*  
321   *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0   *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
322   *     *  
323   */   */
324  void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u)  void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const Paso_Coupler* u_coupler)
325  {  {
326    dim_t n, i, j;    dim_t n, i, j;
327    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
328    double *remote_u;    const double *u=Paso_Coupler_borrowLocalData(u_coupler);
329      const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
330    register double u_i, f_ij;    register double u_i, f_ij;
331    register index_t iptr_ij;    register index_t iptr_ij;
332    
333    n=Paso_SystemMatrix_getTotalNumRows(f);    n=Paso_SystemMatrix_getTotalNumRows(f);
334    pattern=f->pattern;    pattern=f->pattern;
   remote_u=Paso_SystemMatrix_finishCollect(f);  
335    #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)    #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
336    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
337        u_i=u[i];        u_i=u[i];
# Line 429  void Paso_FCTransportProblem_applyPreAnt Line 352  void Paso_FCTransportProblem_applyPreAnt
352    
353    
354  void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,  void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
355                                     const double* QN,const double* QP,double* RN,double* RP)                                     const Paso_Coupler* QN_coupler,const Paso_Coupler* QP_coupler,double* RN,double* RP)
356  {  {
357    dim_t n, i, j;    dim_t n, i, j;
358    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
359      const double *QN=Paso_Coupler_borrowLocalData(QN_coupler);
360      const double *QP=Paso_Coupler_borrowLocalData(QP_coupler);
361    register double f_ij, PP_i, PN_i;    register double f_ij, PP_i, PN_i;
362    register index_t iptr_ij;    register index_t iptr_ij;
363    
# Line 478  void Paso_FCTransportProblem_setRs(const Line 403  void Paso_FCTransportProblem_setRs(const
403    
404  }  }
405    
406  void Paso_FCTransportProblem_addCorrectedFluxes(double* f,Paso_SystemMatrix *flux_matrix,const double* RN,const double* RP)  void Paso_FCTransportProblem_addCorrectedFluxes(double* f,const Paso_SystemMatrix *flux_matrix, const Paso_Coupler* RN_coupler, const Paso_Coupler* RP_coupler)
407  {  {
408    dim_t n, i, j;    dim_t n, i, j;
409    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
   double *remote_RN, *remote_RP;  
410    register double RN_i, RP_i, f_i, f_ij;    register double RN_i, RP_i, f_i, f_ij;
411    register index_t iptr_ij;    register index_t iptr_ij;
412      const double *RN=Paso_Coupler_borrowLocalData(RN_coupler);
413      const double *remote_RN=Paso_Coupler_borrowRemoteData(RN_coupler);
414      const double *RP=Paso_Coupler_borrowLocalData(RP_coupler);
415      const double *remote_RP=Paso_Coupler_borrowRemoteData(RP_coupler);
416    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
417    
   Paso_SystemMatrix_startCollect(flux_matrix,RN);  
418    pattern=flux_matrix->pattern;    pattern=flux_matrix->pattern;
   remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix);  
419    #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)    #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
420    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
421       RN_i=RN[i];       RN_i=RN[i];
# Line 511  void Paso_FCTransportProblem_addCorrecte Line 437  void Paso_FCTransportProblem_addCorrecte
437            f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];            f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];
438            if (f_ij >=0) {            if (f_ij >=0) {
439                f_i+=f_ij*MIN(RP_i,remote_RN[j]);                f_i+=f_ij*MIN(RP_i,remote_RN[j]);
440              }else {
441                  f_i+=f_ij*MIN(RN_i,remote_RP[j]);
442            }            }
443        }        }
444        f[i]+=f_i;        f[i]+=f_i;
445    }    }
   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->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {  
           j=pattern->col_couplePattern->index[iptr_ij];  
           f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];  
           if (f_ij < 0) {  
               f_i+=f_ij*MIN(RN_i,remote_RP[j]);  
           }  
      }  
      f[i]+=f_i;  
   }  
446  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26