/[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 1370 by gross, Wed Jan 2 09:21:43 2008 UTC revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1 
 /* $Id:$ */  
1    
2  /*******************************************************  /*******************************************************
3   *  *
4   *       Copyright 2007 by University of Queensland  * Copyright (c) 2003-2008 by University of Queensland
5   *  * Earth Systems Science Computational Center (ESSCC)
6   *                http://esscc.uq.edu.au  * http://www.uq.edu.au/esscc
7   *        Primary Business: Queensland, Australia  *
8   *  Licensed under the Open Software License version 3.0  * Primary Business: Queensland, Australia
9   *     http://www.opensource.org/licenses/osl-3.0.php  * Licensed under the Open Software License version 3.0
10   *  * http://www.opensource.org/licenses/osl-3.0.php
11   *******************************************************/  *
12    *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 30  Line 31 
31  /**************************************************************/  /**************************************************************/
32    
33  /* free all memory used by                                */  /* free all memory used by                                */
34    
35  void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {  void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {
36       if (in!=NULL) {       if (in!=NULL) {
37          in->reference_counter--;          in->reference_counter--;
38          if (in->reference_counter<=0) {          if (in->reference_counter<=0) {
   
39             Paso_SystemMatrix_free(in->transport_matrix);             Paso_SystemMatrix_free(in->transport_matrix);
40             Paso_SystemMatrix_free(in->flux_matrix);             Paso_SystemMatrix_free(in->mass_matrix);
41               Paso_SystemMatrix_free(in->iteration_matrix);
42             Paso_MPIInfo_free(in->mpi_info);             Paso_MPIInfo_free(in->mpi_info);
43               Paso_Coupler_free(in->u_coupler);
44             MEMFREE(in->u);             MEMFREE(in->u);
            MEMFREE(in->lumped_mass_matrix);  
            MEMFREE(in->row_sum_flux_matrix);  
            MEMFREE(in->transport_matrix_diagonal);  
            MEMFREE(in->colorOf);  
45             MEMFREE(in->main_iptr);             MEMFREE(in->main_iptr);
46               MEMFREE(in->lumped_mass_matrix);
47               MEMFREE(in->main_diagonal_low_order_transport_matrix);
48             MEMFREE(in);             MEMFREE(in);
49          }          }
50      }      }
# Line 60  Paso_FCTransportProblem* Paso_FCTranspor Line 60  Paso_FCTransportProblem* Paso_FCTranspor
60  Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {  Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
61     return in->transport_matrix;     return in->transport_matrix;
62  }  }
63    Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
64  Paso_SystemMatrix* Paso_FCTransportProblem_borrowFluxMatrix(Paso_FCTransportProblem* in) {     return in->mass_matrix;
     return in->flux_matrix;  
65  }  }
66    
67  double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {  double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
# Line 73  dim_t Paso_FCTransportProblem_getTotalNu Line 72  dim_t Paso_FCTransportProblem_getTotalNu
72      return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);      return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
73  }  }
74    
75  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)
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 89  Paso_FCTransportProblem* Paso_FCTranspor Line 86  Paso_FCTransportProblem* Paso_FCTranspor
86    
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         out->reference_counter=0;
90       out->theta=theta;       out->theta=theta;
91       out->dt_max=dt_max;       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;
96         out->u=NULL;
97       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
   
      out->colorOf=NULL;  
98       out->main_iptr=NULL;       out->main_iptr=NULL;
99       out->lumped_mass_matrix=NULL;       out->lumped_mass_matrix=NULL;
100       out->row_sum_flux_matrix=NULL;       out->main_diagonal_low_order_transport_matrix=NULL;
      out->transport_matrix_diagonal=NULL;  
101    
102       if (Paso_noError()) {       if (Paso_noError()) {
103           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
104    
105           out->colorOf=MEMALLOC(n,index_t);           out->u=MEMALLOC(n,double);
106           out->main_iptr=MEMALLOC(n,index_t);           out->main_iptr=MEMALLOC(n,index_t);
107           out->lumped_mass_matrix=MEMALLOC(n,double);           out->lumped_mass_matrix=MEMALLOC(n,double);
108           out->row_sum_flux_matrix=MEMALLOC(n,double);           out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
109           out->transport_matrix_diagonal=MEMALLOC(n,double);           out->u_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(out),block_size);
          out->u=MEMALLOC(n,double);  
110    
111           if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) ||           if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) ||
112                   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)) && Paso_noError()  ) {
113                            
114               printf("Paso_SolverFCT_getFCTransportProblem: Revise coloring!!\n");               #pragma omp parallel
115               Paso_Pattern_color(pattern->mainPattern,&(out->num_colors),out->colorOf);               {
116                     #pragma omp for schedule(static) private(i)
117                                 for (i = 0; i < n; ++i) {
118               #pragma omp parallel for schedule(static) private(i)                      out->lumped_mass_matrix[i]=0.;
119               for (i = 0; i < n; ++i) {                      out->main_diagonal_low_order_transport_matrix[i]=0.;
120                  out->lumped_mass_matrix[i]=0.;                      out->u[i]=0.;
121                  out->row_sum_flux_matrix[i]=0.;                   }
122                  out->u[i]=0.;                   /* identify the main diagonals */
123               }                   #pragma omp for schedule(static) private(i,iptr,iptr_main)
124               /* identify the main diagonals */                   for (i = 0; i < n; ++i) {
125               #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)                          iptr_main=pattern->mainPattern->ptr[0]-1;
126               for (i = 0; i < n; ++i) {                          for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
127                      iptr_main=pattern->mainPattern->ptr[0]-1;                                if (pattern->mainPattern->index[iptr]==i) {
128                      for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {                                     iptr_main=iptr;
129                            if (pattern->mainPattern->index[iptr]==i) {                                     break;
130                                 iptr_main=iptr;                                }
131                                 break;                          }
132                            }                          out->main_iptr[i]=iptr_main;
133                      }                          if (iptr_main==pattern->mainPattern->ptr[0]-1)
134                      out->main_iptr[i]=iptr_main;                               Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
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    if (Paso_noError()) {    if (Paso_noError()) {
142         out->reference_counter=1;
143       return out;       return out;
144    } else {    } else {
145       Paso_FCTransportProblem_free(out);       Paso_FCTransportProblem_free(out);
# Line 155  Paso_FCTransportProblem* Paso_FCTranspor Line 149  Paso_FCTransportProblem* Paso_FCTranspor
149    
150  void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {  void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
151      dim_t i, n;      dim_t i, n;
152          double local_u_min,u_min;
153      n=Paso_FCTransportProblem_getTotalNumRows(in);      n=Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
154      #pragma omp parallel for schedule(static) private(i)      u_min=LARGE_POSITIVE_FLOAT;
155      for (i = 0; i < n; ++i) {      #pragma omp parallel private(local_u_min)
156           in->u[i]=u[i];      {
157             local_u_min=0.;
158             #pragma omp for schedule(static) private(i)
159             for (i=0;i<n;++i) {
160                  in->u[i]=u[i];
161                  local_u_min=MIN(local_u_min,u[i]);
162             }
163             #pragma omp critical
164             {
165                u_min=MIN(u_min,local_u_min);
166             }
167        }
168        #ifdef PASO_MPI
169             local_u_min=u_min;
170             MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
171        #endif
172        if (u_min<0) {
173           Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_checkinSolution: initial guess must be non-negative.");
174      }      }
175  }  }
176    dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
177    {
178       return in->transport_matrix->row_block_size;
179    }
180    
181    Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
182    {
183       return in->transport_matrix->pattern->col_connector;
184    }

Legend:
Removed from v.1370  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26