/[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 1375 by gross, Wed Jan 9 00:15:05 2008 UTC revision 1401 by gross, Fri Jan 25 04:31:18 2008 UTC
# Line 48  Line 48 
48  void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {  void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {
49    dim_t n,i;    dim_t n,i;
50    index_t color, iptr_ij,j,iptr_ji;    index_t color, iptr_ij,j,iptr_ji;
51    register double d_ij;    register double d_ij, sum;
52    
53    if (fc==NULL) return;    if (fc==NULL) return;
54    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
55    
56    #pragma omp parallel private(color)    #pragma omp parallel for private(i,iptr_ij,j,iptr_ji,d_ij,sum)  schedule(static)
57    {    for (i = 0; i < n; ++i) {
58         /* process main block */       sum=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];
59         for (color=0;color<fc->num_colors;++color) {       for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
60             #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)           j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
61             for (i = 0; i < n; ++i) {           if (j!=i) {
62                 if (fc->colorOf[i]==color) {               /* find entry a[j,i] */
63                    fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];               for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
64                    if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
65                    for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {                      d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],
66                       j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];                                            fc->flux_matrix->mainBlock->val[iptr_ji]);
67                       if (j<i) {                      fc->transport_matrix->mainBlock->val[iptr_ij]+=
68                          /* find entry a[j,i] */                                           alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;
69                          for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {                      sum-=d_ij;
70                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {                      break;
71                                  d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);                  }
72  printf("%d %d -> %e\n",i,j,d_ij);               }
73                                  fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;           }
74                                  fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij;       }
75  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]);       /* TODO process couple block */
76                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;  
77                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;       /* update main diagonal */
78                                  break;       fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=sum;
                             }  
                         }  
                      }  
                           
                   }  
                   /* TODO process couple block */  
                }  
            }  
            #pragma omp barrier  
        }  
79    }    }
80    
81  }  }
82    
83  void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {  void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {
84      /*
85       *   sets fa=transport_matrix*u+anti-diffuison_flux(u)
86       */
87    
88    double *remote_u=NULL;    double *remote_u=NULL;
89        
# Line 125  void Paso_FCTransportProblem_setFlux(Pas Line 118  void Paso_FCTransportProblem_setFlux(Pas
118    
119  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {
120    
121    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;    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;
122    index_t color, iptr_ij,j,iptr_ji, i;    index_t color, iptr_ij,j,iptr_ji, i;
123    dim_t n;    dim_t n;
124    
125    
126    if (fc==NULL) return;    if (fc==NULL) return;
127    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
   /* exchange */  
128    
129    
130    #pragma omp parallel private(color)    #pragma omp parallel
131    {    {
132         for (color=0;color<fc->num_colors;++color) {        /*
133             #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)         * calculate the smootness sensors
134             for (i = 0; i < n; ++i) {        */
135                if (fc->colorOf[i]==color) {        #pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,iptr_ij,a_ij,j,d)
136                    u_i=u[i];        for (i = 0; i < n; ++i) {
137                    /* gather the smoothness sensor */            u_i=u[i];
138                    P_p=0.;            P_p=0.;
139                    P_n=0.;            P_n=0.;
140                    Q_p=0.;            Q_p=0.;
141                    Q_n=0.;            Q_n=0.;
142                    r_p=0.;            #pragma ivdep
143                    r_n=0.;        for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
144                    /* #pragma ivdep */                 a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
145                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];
146                        a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];                 d=u[j]-u_i;
147                        j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];                 if (a_ij<0.) {
148                        d=u[j]-u_i;                    if (d<0.) {
149  printf("%d %d : %e %e :: %e \n",i,j,u_i,u[j],a_ij);                        P_p+=a_ij*d;
150                        if (a_ij<0.) {                    } else {
151                           if (d<0.) {                        P_n+=a_ij*d;
152                              P_p+=a_ij*d;                    }
153                           } else {                 } else {
154                              P_n+=a_ij*d;                    if (d>0.) {
155                           }                      Q_p+=a_ij*d;
156                        } else {                    } else {
157                           if (d>0.) {                      Q_n+=a_ij*d;
158                              Q_p+=a_ij*d;                    }
159                           } else {                 }
160                              Q_n+=a_ij*d;        }
161                           }            #pragma ivdep
162                        }        for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
163                }                 a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
164                    #pragma ivdep                 j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
165                for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {                 d=u_remote[j]-u_i;
166                        a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];                 if (a_ij<0.) {
167                        j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];                     if (d<0.) {
168                        d=u_remote[j]-u_i;                       P_p+=a_ij*d;
169                       } else {
170                        if (a_ij<0.) {                       P_n+=a_ij*d;
171                           if (d<0.) {                     }
172                              P_p+=a_ij*d;                 } else {
173                           } else {                    if (d>0.) {
174                              P_n+=a_ij*d;                       Q_p+=a_ij*d;
175                           }                    } else {
176                        } else {                       Q_n+=a_ij*d;
177                           if (d>0.) {                    }
178                              Q_p+=a_ij*d;                 }
179                           } else {        }
180                              Q_n+=a_ij*d;            /* set the smoothness indicators */
181                           }            fc->r_p[i] = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
182              fc->r_n[i] = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.;
183    
184          } /* end of row loop for smootheness indicator */
185    
186          /*
187           * calculate antidiffusion
188          */
189          #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)
190          for (i = 0; i < n; ++i) {
191              u_i=u[i];
192              sum=0;
193              /* anti diffusive flux from main block */
194              for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
195                  a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
196                  j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
197                  if ( i!=j ) {
198                       /* find entry a[j,i] */
199                       for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
200                          if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
201                               a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
202                               d_ij=-MIN3(0,a_ji,a_ij);
203                               if (d_ij > 0.) {
204                                    u_j=u[j];
205                                    if (a_ji >= a_ij) {
206                                        r_ij = u_i>u_j ? fc->r_p[i] : fc->r_n[i];
207                                        f_ij =MIN(r_ij*d_ij,a_ji+d_ij);
208                                    } else {
209                                        r_ij = u_j>u_i ? fc->r_p[j] : fc->r_n[j];
210                                        f_ij =MIN(r_ij*d_ij,a_ij+d_ij);
211                                    }
212                                    sum+=f_ij*(u_i-u_j);
213                               }
214                               break;
215                        }                        }
216                }                     }
217                    /* set the smoothness indicators */                }
218                    r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;            }
219                    r_n = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.;            /* anti diffusive flux from couple block */
220  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);            /* TODO */
                   /* 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]);  
221    
222    
                                     
                                 }  
                                 break;  
                             }  
                         }  
                      }  
                   }  
                   /* anti diffusive flux from couple block */  
223    
224                    /* TODO */            fa[i]+=sum;
225                }        }
226             }    } /* end of parallel block */
            #pragma omp barrier  
        }  
   }  
227  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26