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

temp/paso/src/SolverFCT.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC trunk/paso/src/SolverFCT.c revision 1628 by phornby, Fri Jul 11 13:12:46 2008 UTC
# Line 34  void Paso_FCTransportProblem_free(Paso_F Line 34  void Paso_FCTransportProblem_free(Paso_F
34       if (in!=NULL) {       if (in!=NULL) {
35          in->reference_counter--;          in->reference_counter--;
36          if (in->reference_counter<=0) {          if (in->reference_counter<=0) {
   
37             Paso_SystemMatrix_free(in->transport_matrix);             Paso_SystemMatrix_free(in->transport_matrix);
38             Paso_SystemMatrix_free(in->flux_matrix);             Paso_SystemMatrix_free(in->mass_matrix);
39               Paso_SystemMatrix_free(in->iteration_matrix);
40             Paso_MPIInfo_free(in->mpi_info);             Paso_MPIInfo_free(in->mpi_info);
   
41             MEMFREE(in->u);             MEMFREE(in->u);
42             MEMFREE(in->lumped_mass_matrix);             Paso_Coupler_free(in->u_coupler);
            MEMFREE(in->row_sum_flux_matrix);  
            MEMFREE(in->transport_matrix_diagonal);  
            MEMFREE(in->colorOf);  
43             MEMFREE(in->main_iptr);             MEMFREE(in->main_iptr);
44               MEMFREE(in->lumped_mass_matrix);
45               MEMFREE(in->main_diagonal_low_order_transport_matrix);
46             MEMFREE(in);             MEMFREE(in);
47          }          }
48      }      }
# Line 60  Paso_FCTransportProblem* Paso_FCTranspor Line 58  Paso_FCTransportProblem* Paso_FCTranspor
58  Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {  Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
59     return in->transport_matrix;     return in->transport_matrix;
60  }  }
61    Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
62  Paso_SystemMatrix* Paso_FCTransportProblem_borrowFluxMatrix(Paso_FCTransportProblem* in) {     return in->mass_matrix;
     return in->flux_matrix;  
63  }  }
64    
65  double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {  double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
# Line 73  dim_t Paso_FCTransportProblem_getTotalNu Line 70  dim_t Paso_FCTransportProblem_getTotalNu
70      return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);      return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
71  }  }
72    
73  Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, double dt_max, Paso_SystemMatrixPattern *pattern, int block_size  Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size
74    
75    
76  ) {  ) {
77       Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1;  /* at the moment only block size 1 is supported */       Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1;  /* at the moment only block size 1 is supported */
78       Paso_FCTransportProblem* out=NULL;       Paso_FCTransportProblem* out=NULL;
79       dim_t n,i;       dim_t n,i;
80       index_t iptr,iptr_main,k;       index_t iptr,iptr_main;
81    
82       if ((theta<0.) || (theta >1.)) {       if ((theta<0.) || (theta >1.)) {
83          Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");          Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
# Line 90  Paso_FCTransportProblem* Paso_FCTranspor Line 87  Paso_FCTransportProblem* Paso_FCTranspor
87       out=MEMALLOC(1,Paso_FCTransportProblem);       out=MEMALLOC(1,Paso_FCTransportProblem);
88       if (Paso_checkPtr(out)) return NULL;       if (Paso_checkPtr(out)) return NULL;
89    
90    
91       out->theta=theta;       out->theta=theta;
92       out->dt_max=dt_max;       out->dt_max=LARGE_POSITIVE_FLOAT;
93       out->valid_matrices=FALSE;       out->valid_matrices=FALSE;
94       out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);       out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
95       Paso_SystemMatrix_allocBuffer(out->transport_matrix);       out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
96       out->flux_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);       out->iteration_matrix=NULL;
97         out->u=NULL;
98         out->u_coupler=Paso_Coupler_alloc(pattern->col_connector,block_size);
99         out->u_min=0.;
100       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
   
      out->colorOf=NULL;  
101       out->main_iptr=NULL;       out->main_iptr=NULL;
102       out->lumped_mass_matrix=NULL;       out->lumped_mass_matrix=NULL;
103       out->row_sum_flux_matrix=NULL;       out->main_diagonal_low_order_transport_matrix=NULL;
      out->transport_matrix_diagonal=NULL;  
104    
105       if (Paso_noError()) {       if (Paso_noError()) {
106           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
107    
108           out->colorOf=MEMALLOC(n,index_t);           out->u=MEMALLOC(n,double);
109           out->main_iptr=MEMALLOC(n,index_t);           out->main_iptr=MEMALLOC(n,index_t);
110           out->lumped_mass_matrix=MEMALLOC(n,double);           out->lumped_mass_matrix=MEMALLOC(n,double);
111           out->row_sum_flux_matrix=MEMALLOC(n,double);           out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
          out->transport_matrix_diagonal=MEMALLOC(n,double);  
          out->u=MEMALLOC(n,double);  
112    
113           if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) ||           if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) ||
114                   Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->transport_matrix_diagonal) || Paso_checkPtr(out->row_sum_flux_matrix) || Paso_checkPtr(out->u)) ) {                   Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix))  ) {
115                            
116               printf("Paso_SolverFCT_getFCTransportProblem: Revise coloring!!\n");               #pragma omp parallel
117               Paso_Pattern_color(pattern->mainPattern,&(out->num_colors),out->colorOf);               {
118                     #pragma omp for schedule(static) private(i)
119                                 for (i = 0; i < n; ++i) {
120               #pragma omp parallel for schedule(static) private(i)                      out->lumped_mass_matrix[i]=0.;
121               for (i = 0; i < n; ++i) {                      out->main_diagonal_low_order_transport_matrix[i]=0.;
122                  out->lumped_mass_matrix[i]=0.;                      out->u[i]=0.;
123                  out->row_sum_flux_matrix[i]=0.;                   }
124                  out->u[i]=0.;                   /* identify the main diagonals */
125               }                   #pragma omp for schedule(static) private(i,iptr,iptr_main)
126               /* identify the main diagonals */                   for (i = 0; i < n; ++i) {
127               #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)                          iptr_main=pattern->mainPattern->ptr[0]-1;
128               for (i = 0; i < n; ++i) {                          for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
129                      iptr_main=pattern->mainPattern->ptr[0]-1;                                if (pattern->mainPattern->index[iptr]==i) {
130                      for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {                                     iptr_main=iptr;
131                            if (pattern->mainPattern->index[iptr]==i) {                                     break;
132                                 iptr_main=iptr;                                }
133                                 break;                          }
134                            }                          out->main_iptr[i]=iptr_main;
135                      }                          if (iptr_main==pattern->mainPattern->ptr[0]-1)
136                      out->main_iptr[i]=iptr_main;                               Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
137                      if (iptr_main==pattern->mainPattern->ptr[0]-1)                   }
138                           Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");      
139               }               }
   
140        }        }
141    
142    }    }
# Line 155  Paso_FCTransportProblem* Paso_FCTranspor Line 150  Paso_FCTransportProblem* Paso_FCTranspor
150    
151  void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {  void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
152      dim_t i, n;      dim_t i, n;
153          double local_u_min,u_min;
154      n=Paso_FCTransportProblem_getTotalNumRows(in);      n=Paso_FCTransportProblem_getTotalNumRows(in);
155      
156        u_min=LARGE_POSITIVE_FLOAT;
157        #pragma omp parallel private(local_u_min)
158        {
159             local_u_min=0.;
160             #pragma omp for schedule(static) private(i)
161             for (i=0;i<n;++i) local_u_min=MIN(local_u_min,u[i]);
162             #pragma omp critical
163             {
164                u_min=MIN(u_min,local_u_min);
165             }
166        }
167        #ifdef PASO_MPI
168             local_u_min=u_min;
169             MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
170        #endif
171        in->u_min=u_min*0;
172      #pragma omp parallel for schedule(static) private(i)      #pragma omp parallel for schedule(static) private(i)
173      for (i = 0; i < n; ++i) {      for (i=0;i<n;++i) {
174           in->u[i]=u[i];          in->u[i]=u[i]-(in->u_min);
175      }      }
176        Paso_Coupler_startCollect(in->u_coupler,in->u);
177        Paso_Coupler_finishCollect(in->u_coupler);
178    }
179    dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
180    {
181       return in->transport_matrix->row_block_size;
182    }
183    
184    Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
185    {
186       return in->transport_matrix->pattern->col_connector;
187  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26