/[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 1362 by gross, Mon Dec 17 02:28:16 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);
42             MEMFREE(in->lumped_mass_matrix);             Paso_Coupler_free(in->u_coupler);
            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 51  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 76  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 86  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;
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->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->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
112    
113           if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->lumped_mass_matrix)) ) {           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))  ) {
              printf("Paso_SolverFCT_getFCTransportProblem: Revise coloring!!\n");  
              Paso_Pattern_color(pattern->mainPattern,&(out->num_colors),out->colorOf);  
   
115                            
116               /* identify the main diagonals */               #pragma omp parallel
117               #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)               {
118               for (i = 0; i < n; ++i) {                   #pragma omp for schedule(static) private(i)
119                  for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; ++iptr) {                   for (i = 0; i < n; ++i) {
120                      iptr_main=pattern->mainPattern->ptr[0]-1;                      out->lumped_mass_matrix[i]=0.;
121                      for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {                      out->main_diagonal_low_order_transport_matrix[i]=0.;
122                            if (pattern->mainPattern->index[iptr]==i) {                      out->u[i]=0.;
123                                 iptr_main=iptr;                   }
124                                 break;                   /* identify the main diagonals */
125                            }                   #pragma omp for schedule(static) private(i,iptr,iptr_main)
126                      }                   for (i = 0; i < n; ++i) {
127                      out->main_iptr[i]=iptr_main;                          iptr_main=pattern->mainPattern->ptr[0]-1;
128                      if (iptr_main==pattern->mainPattern->ptr[0]-1)                          for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
129                           Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");                                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    }    }
# Line 136  Paso_FCTransportProblem* Paso_FCTranspor Line 147  Paso_FCTransportProblem* Paso_FCTranspor
147       return NULL;       return NULL;
148    }    }
149  }  }
150    
151    void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
152        dim_t i, n;
153        double local_u_min,u_min;
154        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)
173        for (i=0;i<n;++i) {
174            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.1362  
changed lines
  Added in v.1628

  ViewVC Help
Powered by ViewVC 1.1.26