/[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

trunk/paso/src/Solver_FluxControl.c revision 1361 by gross, Fri Dec 14 09:26:51 2007 UTC trunk/paso/src/SolverFCT_FluxControl.c revision 1370 by gross, Wed Jan 2 09:21:43 2008 UTC
# Line 36  Line 36 
36    
37  /**************************************************************/  /**************************************************************/
38    
 /* free all memory used by FluxControl                                */  
   
 void Paso_Solver_FluxControl_free(Paso_Solver_FluxControl* in) {  
      if (in!=NULL) {  
         Paso_SystemMatrix_freeBuffer(in->matrix);  
         Paso_SystemMatrix_free(in->matrix);  
         MEMFREE(in->colorOf);  
         MEMFREE(in->main_iptr);  
         MEMFREE(in);  
      }  
 }  
   
 /**************************************************************/  
   
 /*   constructs a flux control mechanism                      */  
   
 Paso_Solver_FluxControl* Paso_SolverFCT_getFluxControl(Paso_SystemMatrix * A) {  
   
   Paso_Solver_FluxControl* out=NULL;  
   dim_t n,i;  
   index_t iptr,iptr_main,k;  
   
   if (A==NULL) return out;  
   n=Paso_SystemMatrix_getTotalNumRows(A);  
   if (A->block_size!=1) {  
         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_getFluxControl: block size > 1 is not supported.");  
         return NULL;  
   }  
   out=MEMALLOC(1,Paso_Solver_FluxControl);  
   if (Paso_checkPtr(out)) return NULL;  
   
   out->matrix=Paso_SystemMatrix_reference(A);  
   out->colorOf=NULL;  
   out->main_iptr=NULL;  
     
   
   /* allocations: */    
   out->colorOf=MEMALLOC(n,index_t);  
   out->main_iptr=MEMALLOC(n,index_t);  
   if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) ) ) {  
       printf("Paso_SolverFCT_getFluxControl: Revise coloring!!\n");  
       Paso_Pattern_color(A->mainBlock->pattern,&(out->num_colors),out->colorOf);  
       Paso_SystemMatrix_allocBuffer(A);  
   
       #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)  
       for (i = 0; i < n; ++i) {  
         for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {  
              iptr_main=A->mainBlock->pattern->ptr[0]-1;  
               for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; iptr++) {  
                    if (A->mainBlock->pattern->index[iptr]==i) {  
                         iptr_main=iptr;  
                         break;  
                    }  
                }  
                out->main_iptr[i]=iptr_main;  
                if (iptr_main==A->mainBlock->pattern->ptr[0]-1)  
                   Paso_setError(VALUE_ERROR, "Paso_SolverFCT_getFluxControl: no main diagonal");  
            }  
        }  
   
   }  
   if (Paso_noError()) {  
      return out;  
   } else {  
      Paso_Solver_FluxControl_free(out);  
      return NULL;  
   }  
 }  
   
 /**************************************************************/  
   
39  /* adds A plus stabelising diffusion into the matrix B        */  /* adds A plus stabelising diffusion into the matrix B        */
40    
41  /* d_ij=alpha*max(0,-a[i,j],-a[j,i])  */  /* d_ij=alpha*max(0,-a[i,j],-a[j,i])  */
# Line 115  Paso_Solver_FluxControl* Paso_SolverFCT_ Line 44  Paso_Solver_FluxControl* Paso_SolverFCT_
44  /* b[i,i]-=alpha*d_ij  */  /* b[i,i]-=alpha*d_ij  */
45  /* b[j,j]-=alpha*d_ij  */  /* b[j,j]-=alpha*d_ij  */
46    
47  void Paso_Solver_FluxControl_addDiffusion(Paso_Solver_FluxControl * fc, double alpha, Paso_SystemMatrix * B) {  void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {
48    dim_t n,i;    dim_t n,i;
49    index_t color, iptr_ij,j,iptr_ji;    index_t color, iptr_ij,j,iptr_ji;
50    register double d_ij;    register double d_ij;
51    
52    if (fc==NULL) return;    if (fc==NULL) return;
53    n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
   /* TODO test - same pattern + block size */  
54    
55    #pragma omp parallel private(color)    #pragma omp parallel private(color)
56    {    {
# Line 131  void Paso_Solver_FluxControl_addDiffusio Line 59  void Paso_Solver_FluxControl_addDiffusio
59             #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)             #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)
60             for (i = 0; i < n; ++i) {             for (i = 0; i < n; ++i) {
61                 if (fc->colorOf[i]==color) {                 if (fc->colorOf[i]==color) {
62                    for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {                    fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];
63                       j=fc->matrix->mainBlock->pattern->index[iptr_ij];  
64                       if (i<j) {                    for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
65                         j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
66                         if (j<i) {
67                          /* find entry a[j,i] */                          /* find entry a[j,i] */
68                          for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->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) {
69                              if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
70                                  d_ij=(-alpha)*MIN3(0.,fc->matrix->mainBlock->val[iptr_ij],fc->matrix->mainBlock->val[iptr_ji]);                                  d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);
71                                  B->mainBlock->val[iptr_ij]+=alpha*fc->matrix->mainBlock->val[iptr_ij]+d_ij;  printf("%d %d -> %e\n",i,j,d_ij);
72                                  B->mainBlock->val[iptr_ji]+=alpha*fc->matrix->mainBlock->val[iptr_ji]+d_ij;                                  fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;
73                                  B->mainBlock->val[fc->main_iptr[i]]-=d_ij;                                  fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij;
74                                  B->mainBlock->val[fc->main_iptr[j]]-=d_ij;  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]);
75                                    fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;
76                                    fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;
77                                  break;                                  break;
78                              }                              }
79                          }                          }
# Line 156  void Paso_Solver_FluxControl_addDiffusio Line 88  void Paso_Solver_FluxControl_addDiffusio
88    }    }
89    
90  }  }
91    
92    void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {
93    
94      double *remote_u=NULL;
95      
96      if (fc==NULL) return;
97      Paso_SystemMatrix_startCollect(fc->transport_matrix,u);
98      /* process main block */
99      Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);
100      /* finish exchange */
101      remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);
102      /* process couple block */
103      Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);
104    
105      Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);
106    }
107  /**************************************************************/  /**************************************************************/
108    
109  /* adds antidiffusion to fa    /* adds antidiffusion to fa  
# Line 174  void Paso_Solver_FluxControl_addDiffusio Line 122  void Paso_Solver_FluxControl_addDiffusio
122    
123  */  */
124    
125  void Paso_Solver_FluxControl_setAntiDiffusiveFlux(Paso_Solver_FluxControl * fc, double * u, double* fa) {  void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {
126    
127    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;
   double *u_remote=NULL;  
128    index_t color, iptr_ij,j,iptr_ji, i;    index_t color, iptr_ij,j,iptr_ji, i;
129    dim_t n;    dim_t n;
130    
131    
132    if (fc==NULL) return;    if (fc==NULL) return;
133    n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);    n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
134    /* exchange */    /* exchange */
135    Paso_SystemMatrix_startCollect(fc->matrix,u);  
   u_remote=Paso_SystemMatrix_finishCollect(fc->matrix);  
136    
137    #pragma omp parallel private(color)    #pragma omp parallel private(color)
138    {    {
# Line 201  void Paso_Solver_FluxControl_setAntiDiff Line 147  void Paso_Solver_FluxControl_setAntiDiff
147                    Q_p=0.;                    Q_p=0.;
148                    Q_n=0.;                    Q_n=0.;
149                    #pragma ivdep                    #pragma ivdep
150                for (iptr_ij=(fc->matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->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) {
151                        a_ij=fc->matrix->mainBlock->val[iptr_ij];                        a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
152                        j=fc->matrix->mainBlock->pattern->index[iptr_ij];                        j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
153                        d=u[j]-u_i;                        d=u[j]-u_i;
154                        if (a_ij<0.) {                        if (a_ij<0.) {
155                           if (d<0.) {                           if (d<0.) {
# Line 220  void Paso_Solver_FluxControl_setAntiDiff Line 166  void Paso_Solver_FluxControl_setAntiDiff
166                        }                        }
167                }                }
168                    #pragma ivdep                    #pragma ivdep
169                for (iptr_ij=(fc->matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {                for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
170                        a_ij=fc->matrix->coupleBlock->val[iptr_ij];                        a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
171                        j=fc->matrix->coupleBlock->pattern->index[iptr_ij];                        j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
172                        d=u_remote[j]-u_i;                        d=u_remote[j]-u_i;
173                        if (a_ij<0.) {                        if (a_ij<0.) {
174                           if (d<0.) {                           if (d<0.) {
# Line 242  void Paso_Solver_FluxControl_setAntiDiff Line 188  void Paso_Solver_FluxControl_setAntiDiff
188                    r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;                    r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
189                    r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;                    r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;
190                    /* anti diffusive flux from main block */                    /* anti diffusive flux from main block */
191                    for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->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) {
192                       a_ij=fc->matrix->mainBlock->val[iptr_ij];                       a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
193                       j=fc->matrix->mainBlock->pattern->index[iptr_ij];                       j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
194                       if (a_ij < 0 && i!=j) {                       if (a_ij < 0 && i!=j) {
195                          /* find entry a[j,i] */                          /* find entry a[j,i] */
196                          for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {                          for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
197                              if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {                              if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
198                                  a_ji=fc->matrix->mainBlock->val[iptr_ji];                                  a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
199                                  if  (a_ji > a_ij || (a_ji == a_ij && j<i) ) {                                  if  (a_ji > a_ij || (a_ji == a_ij && j<i) ) {
200                                      u_j=u[j];                                      u_j=u[j];
201                                      r_ij = u_i>u_j ? r_p : r_n;                                      r_ij = u_i>u_j ? r_p : r_n;

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

  ViewVC Help
Powered by ViewVC 1.1.26