/[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 1415 by ksteube, Thu Feb 21 04:57:17 2008 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->mass_matrix);             Paso_SystemMatrix_free(in->mass_matrix);
39             Paso_SystemMatrix_free(in->iteration_matrix);             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               Paso_Coupler_free(in->u_coupler);
43             MEMFREE(in->main_iptr);             MEMFREE(in->main_iptr);
44             MEMFREE(in->lumped_mass_matrix);             MEMFREE(in->lumped_mass_matrix);
45             MEMFREE(in->main_diagonal_low_order_transport_matrix);             MEMFREE(in->main_diagonal_low_order_transport_matrix);
# Line 77  Paso_FCTransportProblem* Paso_FCTranspor Line 77  Paso_FCTransportProblem* Paso_FCTranspor
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 92  Paso_FCTransportProblem* Paso_FCTranspor Line 92  Paso_FCTransportProblem* Paso_FCTranspor
92       out->dt_max=LARGE_POSITIVE_FLOAT;       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);
      Paso_SystemMatrix_allocBuffer(out->transport_matrix);  
95       out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);       out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
96       out->iteration_matrix=NULL;       out->iteration_matrix=NULL;
   
      out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);  
97       out->u=NULL;       out->u=NULL;
98         out->u_coupler=Paso_Coupler_alloc(pattern->col_connector,block_size);
99       out->u_min=0.;       out->u_min=0.;
100         out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
101       out->main_iptr=NULL;       out->main_iptr=NULL;
102       out->lumped_mass_matrix=NULL;       out->lumped_mass_matrix=NULL;
103       out->main_diagonal_low_order_transport_matrix=NULL;       out->main_diagonal_low_order_transport_matrix=NULL;
# Line 114  Paso_FCTransportProblem* Paso_FCTranspor Line 113  Paso_FCTransportProblem* Paso_FCTranspor
113           if ( ! (Paso_checkPtr(out->u) || 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->main_diagonal_low_order_transport_matrix))  ) {                   Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix))  ) {
115                            
116               #pragma omp parallel for schedule(static) private(i)               #pragma omp parallel
117               for (i = 0; i < n; ++i) {               {
118                  out->lumped_mass_matrix[i]=0.;                   #pragma omp for schedule(static) private(i)
119                  out->main_diagonal_low_order_transport_matrix[i]=0.;                   for (i = 0; i < n; ++i) {
120                  out->u[i]=0.;                      out->lumped_mass_matrix[i]=0.;
121               }                      out->main_diagonal_low_order_transport_matrix[i]=0.;
122               /* identify the main diagonals */                      out->u[i]=0.;
123               #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)                   }
124               for (i = 0; i < n; ++i) {                   /* identify the main diagonals */
125                      iptr_main=pattern->mainPattern->ptr[0]-1;                   #pragma omp for schedule(static) private(i,iptr,iptr_main)
126                      for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {                   for (i = 0; i < n; ++i) {
127                            if (pattern->mainPattern->index[iptr]==i) {                          iptr_main=pattern->mainPattern->ptr[0]-1;
128                                 iptr_main=iptr;                          for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
129                                 break;                                if (pattern->mainPattern->index[iptr]==i) {
130                            }                                     iptr_main=iptr;
131                      }                                     break;
132                      out->main_iptr[i]=iptr_main;                                }
133                      if (iptr_main==pattern->mainPattern->ptr[0]-1)                          }
134                           Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");                          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 149  Paso_FCTransportProblem* Paso_FCTranspor Line 151  Paso_FCTransportProblem* Paso_FCTranspor
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;      double local_u_min,u_min;
     
154      n=Paso_FCTransportProblem_getTotalNumRows(in);      n=Paso_FCTransportProblem_getTotalNumRows(in);
155      
156      u_min=LARGE_POSITIVE_FLOAT;      u_min=LARGE_POSITIVE_FLOAT;
157      #pragma omp parallel private(local_u_min)      #pragma omp parallel private(local_u_min)
158      {      {
# Line 166  void Paso_FCTransportProblem_checkinSolu Line 168  void Paso_FCTransportProblem_checkinSolu
168           local_u_min=u_min;           local_u_min=u_min;
169           MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);           MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
170      #endif      #endif
171      in->u_min=u_min;      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]-u_min;          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.1415  
changed lines
  Added in v.1628

  ViewVC Help
Powered by ViewVC 1.1.26