/[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 1400 by gross, Thu Jan 24 06:04:31 2008 UTC
# Line 29  Line 29 
29  #define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.)  #define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.)
30  #define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b)))  #define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b)))
31  #define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b))))  #define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b))))
32    #define MC(a,b) (FLUX_S(a,b)*MIN3(ABS((a)+(b))/2,2*ABS(a),2*ABS(b)))
33    
34  #define FLUX_L(a,b) SUPERBEE(a,b)  /* alter for other flux limiter */  #define FLUX_L(a,b) MC(a,b)  /* alter for other flux limiter */
35    
36  #define FLUX_LIMITER(a) FLUX_L(a,1)  #define FLUX_LIMITER(a) FLUX_L(1,a)
37    
38  /**************************************************************/  /**************************************************************/
39    
# Line 68  void Paso_FCTransportProblem_addAdvectiv Line 69  void Paso_FCTransportProblem_addAdvectiv
69                          for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {                          for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
70                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
71                                  d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);                                  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);  /* 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;                                  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;                                  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]);  /* 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]); */
76                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;
77                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;                                  fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;
78                                  break;                                  break;
# Line 90  printf("%d %d -> %e -> %e %e \n",i,j,d_i Line 91  printf("%d %d -> %e -> %e %e \n",i,j,d_i
91  }  }
92    
93  void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {  void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {
94      /*
95       *   sets fa=transport_matrix*u+anti-diffuison_flux(u)
96       */
97    
98    double *remote_u=NULL;    double *remote_u=NULL;
99        
# Line 124  void Paso_FCTransportProblem_setFlux(Pas Line 128  void Paso_FCTransportProblem_setFlux(Pas
128    
129  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) {
130    
131    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;    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;
132    index_t color, iptr_ij,j,iptr_ji, i;    index_t color, iptr_ij,j,iptr_ji, i;
133    dim_t n;    dim_t n;
134    
# Line 137  void Paso_FCTransportProblem_setAntiDiff Line 141  void Paso_FCTransportProblem_setAntiDiff
141    #pragma omp parallel private(color)    #pragma omp parallel private(color)
142    {    {
143         for (color=0;color<fc->num_colors;++color) {         for (color=0;color<fc->num_colors;++color) {
144             #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)             #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)
145             for (i = 0; i < n; ++i) {             for (i = 0; i < n; ++i) {
146                if (fc->colorOf[i]==color) {                if (fc->colorOf[i]==color) {
147                    u_i=u[i];                    u_i=u[i];
# Line 146  void Paso_FCTransportProblem_setAntiDiff Line 150  void Paso_FCTransportProblem_setAntiDiff
150                    P_n=0.;                    P_n=0.;
151                    Q_p=0.;                    Q_p=0.;
152                    Q_n=0.;                    Q_n=0.;
153                    #pragma ivdep                    r_p=0.;
154                      r_n=0.;
155                      /* #pragma ivdep */
156                for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {                for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
157                        a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];                        a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
158                        j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];                        j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
159                        d=u[j]-u_i;                        d=u[j]-u_i;
160    /* printf("%d %d : %e %e :: %e \n",i,j,u_i,u[j],a_ij); */
161                        if (a_ij<0.) {                        if (a_ij<0.) {
162                           if (d<0.) {                           if (d<0.) {
163                              P_p+=a_ij*d;                              P_p+=a_ij*d;
# Line 170  void Paso_FCTransportProblem_setAntiDiff Line 177  void Paso_FCTransportProblem_setAntiDiff
177                        a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];                        a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
178                        j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];                        j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
179                        d=u_remote[j]-u_i;                        d=u_remote[j]-u_i;
180    
181                        if (a_ij<0.) {                        if (a_ij<0.) {
182                           if (d<0.) {                           if (d<0.) {
183                              P_p+=a_ij*d;                              P_p+=a_ij*d;
# Line 186  void Paso_FCTransportProblem_setAntiDiff Line 194  void Paso_FCTransportProblem_setAntiDiff
194                }                }
195                    /* set the smoothness indicators */                    /* set the smoothness indicators */
196                    r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;                    r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
197                    r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;                    r_n = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.;
198    /* printf("Flux control %d: %e %e %e : %e %e %e : %e\n",i,Q_p,P_p,r_p,Q_n,P_n,r_n,u_i);  */
199                    /* anti diffusive flux from main block */                    /* anti diffusive flux from main block */
200                    for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {                    for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
201                       a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];                       a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
202                       j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];                       j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
203                       if (a_ij < 0 && i!=j) {                       if ( i!=j ) {
204                          /* find entry a[j,i] */                          /* find entry a[j,i] */
205                          for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {                          for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
206                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
207                                  a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];                                  a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
208                                  if  (a_ji > a_ij || (a_ji == a_ij && j<i) ) {                                  d_ij=-MIN3(0,a_ji,a_ij);
209                                    if ( (d_ij > 0.) && ( (a_ji > a_ij) || ( (a_ji == a_ij) && (j<i) ) )  ) {
210                                      u_j=u[j];                                      u_j=u[j];
211                                      r_ij = u_i>u_j ? r_p : r_n;                                      r_ij = u_i>u_j ? r_p : r_n;
212                                      f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j);                                      f_ij =MIN(r_ij*d_ij,a_ji+d_ij)*(u_i-u_j);
213    
214                                      fa[i]+=f_ij;                                      fa[i]+=f_ij;
215                                      fa[j]-=f_ij;                                      fa[j]-=f_ij;
216                                      break;  /* 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]); */
217    
218    
219                                      
220                                  }                                  }
221                                    break;
222                              }                              }
223                          }                          }
224                       }                       }

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

  ViewVC Help
Powered by ViewVC 1.1.26