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

revision 1366 by gross, Tue Dec 18 05:49:17 2007 UTC 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->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 53  Paso_FCTransportProblem* Paso_FCTranspor Line 52  Paso_FCTransportProblem* Paso_FCTranspor
52       if (in!=NULL) {       if (in!=NULL) {
53          ++(in->reference_counter);          ++(in->reference_counter);
54       }       }
55         return in;
56  }      }    
57    
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 71  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 88  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=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;
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->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->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               /* identify the main diagonals */                      out->lumped_mass_matrix[i]=0.;
121               #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)                      out->main_diagonal_low_order_transport_matrix[i]=0.;
122               for (i = 0; i < n; ++i) {                      out->u[i]=0.;
123                  for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; ++iptr) {                   }
124                      iptr_main=pattern->mainPattern->ptr[0]-1;                   /* identify the main diagonals */
125                      for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {                   #pragma omp for schedule(static) private(i,iptr,iptr_main)
126                            if (pattern->mainPattern->index[iptr]==i) {                   for (i = 0; i < n; ++i) {
127                                 iptr_main=iptr;                          iptr_main=pattern->mainPattern->ptr[0]-1;
128                                 break;                          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                      out->main_iptr[i]=iptr_main;                                     break;
132                      if (iptr_main==pattern->mainPattern->ptr[0]-1)                                }
133                           Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");                          }
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    }    }
# Line 146  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.1366  
changed lines
  Added in v.1628

  ViewVC Help
Powered by ViewVC 1.1.26