/[escript]/trunk/paso/src/SolverFCT.c
ViewVC logotype

Diff of /trunk/paso/src/SolverFCT.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.c revision 1628 by phornby, Fri Jul 11 13:12:46 2008 UTC
# Line 13  Line 13 
13    
14  /**************************************************************/  /**************************************************************/
15    
16  /* Paso: FluxControl                                          */  /* Paso: FCTransportProblem                                          */
17    
18  /**************************************************************/  /**************************************************************/
19    
# Line 26  Line 26 
26  #include "SolverFCT.h"  #include "SolverFCT.h"
27  #include "PasoUtil.h"  #include "PasoUtil.h"
28    
 #define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.)  
 #define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b)))  
 #define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b))))  
   
 #define FLUX_L(a,b) SUPERBEE(a,b)  /* alter for other flux limiter */  
   
 #define FLUX_LIMITER(a) FLUX_L(a,1)  
29    
30  /**************************************************************/  /**************************************************************/
31    
32  /* free all memory used by FluxControl                                */  /* free all memory used by                                */
33    void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {
34         if (in!=NULL) {
35            in->reference_counter--;
36            if (in->reference_counter<=0) {
37               Paso_SystemMatrix_free(in->transport_matrix);
38               Paso_SystemMatrix_free(in->mass_matrix);
39               Paso_SystemMatrix_free(in->iteration_matrix);
40               Paso_MPIInfo_free(in->mpi_info);
41               MEMFREE(in->u);
42               Paso_Coupler_free(in->u_coupler);
43               MEMFREE(in->main_iptr);
44               MEMFREE(in->lumped_mass_matrix);
45               MEMFREE(in->main_diagonal_low_order_transport_matrix);
46               MEMFREE(in);
47            }
48        }
49    }
50    
51  void Paso_Solver_FluxControl_free(Paso_Solver_FluxControl* in) {  Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
52       if (in!=NULL) {       if (in!=NULL) {
53          Paso_SystemMatrix_freeBuffer(in->matrix);          ++(in->reference_counter);
         Paso_SystemMatrix_free(in->matrix);  
         MEMFREE(in->colorOf);  
         MEMFREE(in->main_iptr);  
         MEMFREE(in);  
54       }       }
55         return in;
56    }    
57    
58    Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
59       return in->transport_matrix;
60    }
61    Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
62       return in->mass_matrix;
63  }  }
64    
65  /**************************************************************/  double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
66        return in->lumped_mass_matrix;
67    }
68    
69    dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
70        return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
71    }
72    
73  /*   constructs a flux control mechanism                      */  Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size
74    
 Paso_Solver_FluxControl* Paso_SolverFCT_getFluxControl(Paso_SystemMatrix * A) {  
75    
76    Paso_Solver_FluxControl* out=NULL;  ) {
77    dim_t n,i;       Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1;  /* at the moment only block size 1 is supported */
78    index_t iptr,iptr_main,k;       Paso_FCTransportProblem* out=NULL;
79         dim_t n,i;
80    if (A==NULL) return out;       index_t iptr,iptr_main;
81    n=Paso_SystemMatrix_getTotalNumRows(A);  
82    if (A->block_size!=1) {       if ((theta<0.) || (theta >1.)) {
83          Paso_setError(TYPE_ERROR,"Paso_SolverFCT_getFluxControl: block size > 1 is not supported.");          Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
84          return NULL;          return NULL;
85    }       }
   out=MEMALLOC(1,Paso_Solver_FluxControl);  
   if (Paso_checkPtr(out)) return NULL;  
86    
87    out->matrix=Paso_SystemMatrix_reference(A);       out=MEMALLOC(1,Paso_FCTransportProblem);
88    out->colorOf=NULL;       if (Paso_checkPtr(out)) return NULL;
89    out->main_iptr=NULL;  
90      
91         out->theta=theta;
92    /* allocations: */         out->dt_max=LARGE_POSITIVE_FLOAT;
93    out->colorOf=MEMALLOC(n,index_t);       out->valid_matrices=FALSE;
94    out->main_iptr=MEMALLOC(n,index_t);       out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
95    if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) ) ) {       out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
96        printf("Paso_SolverFCT_getFluxControl: Revise coloring!!\n");       out->iteration_matrix=NULL;
97        Paso_Pattern_color(A->mainBlock->pattern,&(out->num_colors),out->colorOf);       out->u=NULL;
98        Paso_SystemMatrix_allocBuffer(A);       out->u_coupler=Paso_Coupler_alloc(pattern->col_connector,block_size);
99         out->u_min=0.;
100        #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
101        for (i = 0; i < n; ++i) {       out->main_iptr=NULL;
102          for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {       out->lumped_mass_matrix=NULL;
103               iptr_main=A->mainBlock->pattern->ptr[0]-1;       out->main_diagonal_low_order_transport_matrix=NULL;
104                for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; iptr++) {  
105                     if (A->mainBlock->pattern->index[iptr]==i) {       if (Paso_noError()) {
106                          iptr_main=iptr;           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
107                          break;  
108                     }           out->u=MEMALLOC(n,double);
109                 }           out->main_iptr=MEMALLOC(n,index_t);
110                 out->main_iptr[i]=iptr_main;           out->lumped_mass_matrix=MEMALLOC(n,double);
111                 if (iptr_main==A->mainBlock->pattern->ptr[0]-1)           out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
112                    Paso_setError(VALUE_ERROR, "Paso_SolverFCT_getFluxControl: no main diagonal");  
113             }           if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) ||
114         }                   Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix))  ) {
115                
116                 #pragma omp parallel
117                 {
118                     #pragma omp for schedule(static) private(i)
119                     for (i = 0; i < n; ++i) {
120                        out->lumped_mass_matrix[i]=0.;
121                        out->main_diagonal_low_order_transport_matrix[i]=0.;
122                        out->u[i]=0.;
123                     }
124                     /* identify the main diagonals */
125                     #pragma omp for schedule(static) private(i,iptr,iptr_main)
126                     for (i = 0; i < n; ++i) {
127                            iptr_main=pattern->mainPattern->ptr[0]-1;
128                            for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
129                                  if (pattern->mainPattern->index[iptr]==i) {
130                                       iptr_main=iptr;
131                                       break;
132                                  }
133                            }
134                            out->main_iptr[i]=iptr_main;
135                            if (iptr_main==pattern->mainPattern->ptr[0]-1)
136                                 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
137                     }
138        
139                 }
140          }
141    
142    }    }
143    if (Paso_noError()) {    if (Paso_noError()) {
144       return out;       return out;
145    } else {    } else {
146       Paso_Solver_FluxControl_free(out);       Paso_FCTransportProblem_free(out);
147       return NULL;       return NULL;
148    }    }
149  }  }
150    
151  /**************************************************************/  void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
152        dim_t i, n;
153  /* adds A plus stabelising diffusion into the matrix B        */      double local_u_min,u_min;
154        n=Paso_FCTransportProblem_getTotalNumRows(in);
155  /* d_ij=alpha*max(0,-a[i,j],-a[j,i])  */    
156  /* b[i,j]+=alpha*(a[i,j]+d_ij)  */      u_min=LARGE_POSITIVE_FLOAT;
157  /* b[j,i]+=alpha*(a[j,i]+d_ij)  */      #pragma omp parallel private(local_u_min)
158  /* b[i,i]-=alpha*d_ij  */      {
159  /* b[j,j]-=alpha*d_ij  */           local_u_min=0.;
160             #pragma omp for schedule(static) private(i)
161  void Paso_Solver_FluxControl_addDiffusion(Paso_Solver_FluxControl * fc, double alpha, Paso_SystemMatrix * B) {           for (i=0;i<n;++i) local_u_min=MIN(local_u_min,u[i]);
162    dim_t n,i;           #pragma omp critical
163    index_t color, iptr_ij,j,iptr_ji;           {
164    register double d_ij;              u_min=MIN(u_min,local_u_min);
165             }
166    if (fc==NULL) return;      }
167    n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);      #ifdef PASO_MPI
168    /* TODO test - same pattern + block size */           local_u_min=u_min;
169             MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
170    #pragma omp parallel private(color)      #endif
171    {      in->u_min=u_min*0;
172         /* process main block */      #pragma omp parallel for schedule(static) private(i)
173         for (color=0;color<fc->num_colors;++color) {      for (i=0;i<n;++i) {
174             #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij)  schedule(static)          in->u[i]=u[i]-(in->u_min);
175             for (i = 0; i < n; ++i) {      }
176                 if (fc->colorOf[i]==color) {      Paso_Coupler_startCollect(in->u_coupler,in->u);
177                    for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {      Paso_Coupler_finishCollect(in->u_coupler);
178                       j=fc->matrix->mainBlock->pattern->index[iptr_ij];  }
179                       if (i<j) {  dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
180                          /* find entry a[j,i] */  {
181                          for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {     return in->transport_matrix->row_block_size;
                             if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {  
                                 d_ij=(-alpha)*MIN3(0.,fc->matrix->mainBlock->val[iptr_ij],fc->matrix->mainBlock->val[iptr_ji]);  
                                 B->mainBlock->val[iptr_ij]+=alpha*fc->matrix->mainBlock->val[iptr_ij]+d_ij;  
                                 B->mainBlock->val[iptr_ji]+=alpha*fc->matrix->mainBlock->val[iptr_ji]+d_ij;  
                                 B->mainBlock->val[fc->main_iptr[i]]-=d_ij;  
                                 B->mainBlock->val[fc->main_iptr[j]]-=d_ij;  
                                 break;  
                             }  
                         }  
                      }  
                           
                   }  
                   /* TODO process couple block */  
                }  
            }  
            #pragma omp barrier  
        }  
   }  
   
182  }  }
 /**************************************************************/  
183    
184  /* adds antidiffusion to fa    Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
185    {
186     P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])       return in->transport_matrix->pattern->col_connector;
    P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])    
    Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])    
    Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])    
    d_ij=max(0,-a[i,j],-a[j,i])  
    l_ji=max(0,a[j,i],a[j,i]-a[i,j])  
    if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)  
       r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]  
       f_ij =min(FLUX_LIMITER(r_ij)*d_ij,l_ji) (u[i]-u[j])=min(FLUX_LIMITER(r_ij)*a[i,j],a[j,i]-a[i,j]) (u[i]-u[j])  
       fa[i]+=f_ij  
       fa[j]-=f_ij  
   
 */  
   
 void Paso_Solver_FluxControl_setAntiDiffusiveFlux(Paso_Solver_FluxControl * fc, double * u, double* fa) {  
   
   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;  
   index_t color, iptr_ij,j,iptr_ji, i;  
   dim_t n;  
   
   
   if (fc==NULL) return;  
   n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);  
   /* exchange */  
   Paso_SystemMatrix_startCollect(fc->matrix,u);  
   u_remote=Paso_SystemMatrix_finishCollect(fc->matrix);  
   
   #pragma omp parallel private(color)  
   {  
        for (color=0;color<fc->num_colors;++color) {  
            #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)  
            for (i = 0; i < n; ++i) {  
               if (fc->colorOf[i]==color) {  
                   u_i=u[i];  
                   /* gather the smoothness sensor */  
                   P_p=0.;  
                   P_n=0.;  
                   Q_p=0.;  
                   Q_n=0.;  
                   #pragma ivdep  
               for (iptr_ij=(fc->matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {  
                       a_ij=fc->matrix->mainBlock->val[iptr_ij];  
                       j=fc->matrix->mainBlock->pattern->index[iptr_ij];  
                       d=u[j]-u_i;  
                       if (a_ij<0.) {  
                          if (d<0.) {  
                             P_p+=a_ij*d;  
                          } else {  
                             P_n+=a_ij*d;  
                          }  
                       } else {  
                          if (d>0.) {  
                             Q_p+=a_ij*d;  
                          } else {  
                             Q_n+=a_ij*d;  
                          }  
                       }  
               }  
                   #pragma ivdep  
               for (iptr_ij=(fc->matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {  
                       a_ij=fc->matrix->coupleBlock->val[iptr_ij];  
                       j=fc->matrix->coupleBlock->pattern->index[iptr_ij];  
                       d=u_remote[j]-u_i;  
                       if (a_ij<0.) {  
                          if (d<0.) {  
                             P_p+=a_ij*d;  
                          } else {  
                             P_n+=a_ij*d;  
                          }  
                       } else {  
                          if (d>0.) {  
                             Q_p+=a_ij*d;  
                          } else {  
                             Q_n+=a_ij*d;  
                          }  
                       }  
               }  
                   /* set the smoothness indicators */  
                   r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;  
                   r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;  
                   /* anti diffusive flux from main block */  
                   for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {  
                      a_ij=fc->matrix->mainBlock->val[iptr_ij];  
                      j=fc->matrix->mainBlock->pattern->index[iptr_ij];  
                      if (a_ij < 0 && i!=j) {  
                         /* find entry a[j,i] */  
                         for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {  
                             if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {  
                                 a_ji=fc->matrix->mainBlock->val[iptr_ji];  
                                 if  (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*a_ij,a_ji-a_ij)*(u_i-u_j);  
                                     fa[i]+=f_ij;  
                                     fa[j]-=f_ij;  
                                     break;  
                                 }  
                             }  
                         }  
                      }  
                   }  
                   /* anti diffusive flux from couple block */  
   
                   /* TODO */  
               }  
            }  
            #pragma omp barrier  
        }  
   }  
187  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26