/[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 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC revision 1400 by gross, Thu Jan 24 06:04:31 2008 UTC
# Line 69  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 91  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 154  void Paso_FCTransportProblem_setAntiDiff Line 157  void Paso_FCTransportProblem_setAntiDiff
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);  /* 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 192  printf("%d %d : %e %e :: %e \n",i,j,u_i, Line 195  printf("%d %d : %e %e :: %e \n",i,j,u_i,
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("%d: %e %e %e : %e %e %e : %e\n",i,Q_p,P_p,r_p,Q_n,P_n,r_n,u_i);  /* 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];
# Line 210  printf("%d: %e %e %e : %e %e %e : %e\n", Line 213  printf("%d: %e %e %e : %e %e %e : %e\n",
213    
214                                      fa[i]+=f_ij;                                      fa[i]+=f_ij;
215                                      fa[j]-=f_ij;                                      fa[j]-=f_ij;
216  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]);  /* 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                                                                        

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

  ViewVC Help
Powered by ViewVC 1.1.26