/[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 1363 by gross, Mon Dec 17 03:42:03 2007 UTC revision 1552 by gross, Thu May 8 08:52:41 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->lumped_mass_matrix);  
            MEMFREE(in->row_sum_flux_matrix);  
            MEMFREE(in->colorOf);  
42             MEMFREE(in->main_iptr);             MEMFREE(in->main_iptr);
43               MEMFREE(in->lumped_mass_matrix);
44               MEMFREE(in->main_diagonal_low_order_transport_matrix);
45             MEMFREE(in);             MEMFREE(in);
46          }          }
47      }      }
# Line 52  Paso_FCTransportProblem* Paso_FCTranspor Line 51  Paso_FCTransportProblem* Paso_FCTranspor
51       if (in!=NULL) {       if (in!=NULL) {
52          ++(in->reference_counter);          ++(in->reference_counter);
53       }       }
54         return in;
55  }      }    
56    
57  Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {  Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
58     return in->transport_matrix;     return in->transport_matrix;
59  }  }
60    Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
61  Paso_SystemMatrix* Paso_FCTransportProblem_borrowFluxMatrix(Paso_FCTransportProblem* in) {     return in->mass_matrix;
     return in->flux_matrix;  
62  }  }
63    
64  double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {  double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
# Line 87  Paso_FCTransportProblem* Paso_FCTranspor Line 86  Paso_FCTransportProblem* Paso_FCTranspor
86       out=MEMALLOC(1,Paso_FCTransportProblem);       out=MEMALLOC(1,Paso_FCTransportProblem);
87       if (Paso_checkPtr(out)) return NULL;       if (Paso_checkPtr(out)) return NULL;
88    
89    
90       out->theta=theta;       out->theta=theta;
91         out->dt_max=LARGE_POSITIVE_FLOAT;
92       out->valid_matrices=FALSE;       out->valid_matrices=FALSE;
93       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);
94       Paso_SystemMatrix_allocBuffer(out->transport_matrix);       out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
95       out->flux_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);       out->iteration_matrix=NULL;
      out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);  
96    
97       out->colorOf=NULL;       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
98         out->u=NULL;
99         out->u_min=0.;
100       out->main_iptr=NULL;       out->main_iptr=NULL;
101       out->lumped_mass_matrix=NULL;       out->lumped_mass_matrix=NULL;
102       out->row_sum_flux_matrix=NULL;       out->main_diagonal_low_order_transport_matrix=NULL;
103    
104       if (Paso_noError()) {       if (Paso_noError()) {
105           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
106    
107           out->colorOf=MEMALLOC(n,index_t);           out->u=MEMALLOC(n,double);
108           out->main_iptr=MEMALLOC(n,index_t);           out->main_iptr=MEMALLOC(n,index_t);
109           out->lumped_mass_matrix=MEMALLOC(n,double);           out->lumped_mass_matrix=MEMALLOC(n,double);
110           out->row_sum_flux_matrix=MEMALLOC(n,double);           out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
   
          if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) ||  
                  Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->row_sum_flux_matrix)) ) {  
               
              printf("Paso_SolverFCT_getFCTransportProblem: Revise coloring!!\n");  
              Paso_Pattern_color(pattern->mainPattern,&(out->num_colors),out->colorOf);  
111    
112             if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) ||
113                     Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix))  ) {
114                            
115                 #pragma omp parallel for schedule(static) private(i)
116                 for (i = 0; i < n; ++i) {
117                    out->lumped_mass_matrix[i]=0.;
118                    out->main_diagonal_low_order_transport_matrix[i]=0.;
119                    out->u[i]=0.;
120                 }
121               /* identify the main diagonals */               /* identify the main diagonals */
122               #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)               #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
123               for (i = 0; i < n; ++i) {               for (i = 0; i < n; ++i) {
                 for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; ++iptr) {  
124                      iptr_main=pattern->mainPattern->ptr[0]-1;                      iptr_main=pattern->mainPattern->ptr[0]-1;
125                      for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {                      for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
126                            if (pattern->mainPattern->index[iptr]==i) {                            if (pattern->mainPattern->index[iptr]==i) {
# Line 128  Paso_FCTransportProblem* Paso_FCTranspor Line 131  Paso_FCTransportProblem* Paso_FCTranspor
131                      out->main_iptr[i]=iptr_main;                      out->main_iptr[i]=iptr_main;
132                      if (iptr_main==pattern->mainPattern->ptr[0]-1)                      if (iptr_main==pattern->mainPattern->ptr[0]-1)
133                           Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");                           Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
                 }  
134               }               }
135    
136        }        }
# Line 141  Paso_FCTransportProblem* Paso_FCTranspor Line 143  Paso_FCTransportProblem* Paso_FCTranspor
143       return NULL;       return NULL;
144    }    }
145  }  }
146    
147    void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
148        dim_t i, n;
149        double local_u_min,u_min;
150      
151        n=Paso_FCTransportProblem_getTotalNumRows(in);
152        u_min=LARGE_POSITIVE_FLOAT;
153        #pragma omp parallel private(local_u_min)
154        {
155             local_u_min=0.;
156             #pragma omp for schedule(static) private(i)
157             for (i=0;i<n;++i) local_u_min=MIN(local_u_min,u[i]);
158             #pragma omp critical
159             {
160                u_min=MIN(u_min,local_u_min);
161             }
162        }
163        #ifdef PASO_MPI
164             local_u_min=u_min;
165             MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
166        #endif
167        in->u_min=u_min;
168        #pragma omp parallel for schedule(static) private(i)
169        for (i=0;i<n;++i) {
170            in->u[i]=u[i]-u_min;
171        }
172    }

Legend:
Removed from v.1363  
changed lines
  Added in v.1552

  ViewVC Help
Powered by ViewVC 1.1.26