/[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 2196 by gross, Wed Oct 8 03:03:37 2008 UTC revision 2197 by gross, Thu Jan 8 05:49:16 2009 UTC
# Line 42  void Paso_FCTransportProblem_free(Paso_F Line 42  void Paso_FCTransportProblem_free(Paso_F
42             Paso_MPIInfo_free(in->mpi_info);             Paso_MPIInfo_free(in->mpi_info);
43             Paso_Coupler_free(in->u_coupler);             Paso_Coupler_free(in->u_coupler);
44             MEMFREE(in->u);             MEMFREE(in->u);
45               MEMFREE(in->constraint_weights);
46             MEMFREE(in->main_iptr);             MEMFREE(in->main_iptr);
47             MEMFREE(in->lumped_mass_matrix);             MEMFREE(in->lumped_mass_matrix);
48             MEMFREE(in->main_diagonal_low_order_transport_matrix);             MEMFREE(in->main_diagonal_low_order_transport_matrix);
# Line 90  Paso_FCTransportProblem* Paso_FCTranspor Line 91  Paso_FCTransportProblem* Paso_FCTranspor
91       out->reference_counter=0;       out->reference_counter=0;
92       out->theta=theta;       out->theta=theta;
93       out->dt_max=LARGE_POSITIVE_FLOAT;       out->dt_max=LARGE_POSITIVE_FLOAT;
94         out->constraint_factor=sqrt(LARGE_POSITIVE_FLOAT);
95         if (out->theta < 1.) {
96                out->dt_factor=MIN(1./(1.-out->theta),DT_FACTOR_MAX);
97         } else {
98                out->dt_factor=DT_FACTOR_MAX;
99         }
100       out->valid_matrices=FALSE;       out->valid_matrices=FALSE;
101       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);
102       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);
103       out->iteration_matrix=NULL;       out->iteration_matrix=NULL;
104       out->u=NULL;       out->u=NULL;
105         out->constraint_weights=NULL;
106       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
107       out->main_iptr=NULL;       out->main_iptr=NULL;
108       out->lumped_mass_matrix=NULL;       out->lumped_mass_matrix=NULL;
# Line 103  Paso_FCTransportProblem* Paso_FCTranspor Line 111  Paso_FCTransportProblem* Paso_FCTranspor
111       if (Paso_noError()) {       if (Paso_noError()) {
112           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);           n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
113           transport_pattern=out->transport_matrix->pattern;           transport_pattern=out->transport_matrix->pattern;
   
114           out->u=MEMALLOC(n,double);           out->u=MEMALLOC(n,double);
115             out->constraint_weights=MEMALLOC(n,double);
116           out->main_iptr=MEMALLOC(n,index_t);           out->main_iptr=MEMALLOC(n,index_t);
117           out->lumped_mass_matrix=MEMALLOC(n,double);           out->lumped_mass_matrix=MEMALLOC(n,double);
118           out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);           out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
119           out->u_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(out),block_size);           out->u_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(out),block_size);
120    
121           if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) ||           if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->constraint_weights) ||
122                   Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Paso_noError()  ) {                   Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Paso_noError()  ) {
123                            
124               #pragma omp parallel               #pragma omp parallel
# Line 190  index_t Paso_FCTransportProblem_getTypeI Line 198  index_t Paso_FCTransportProblem_getTypeI
198     return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;     return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
199  }  }
200    
201    void Paso_FCTransportProblem_setUpConstraint(Paso_FCTransportProblem* fctp,  const double* q, const double factor)
202    {
203       dim_t i, n;
204       register double m, rtmp;
205       double factor2= fctp->dt_factor * factor;
206       n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
207    
208       if ( fctp->valid_matrices ) {
209          Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_insertConstraint: you must not insert a constraint is a valid system.");
210          return;
211       }
212       if (factor<=0) {
213          Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_insertConstraint: constraint_factor needs to be positive.");
214          return;
215       }
216    
217      
218       #pragma omp for schedule(static) private(i,m,rtmp)
219       for (i=0;i<n;++i) {
220            if (q[i]>0) {
221               m=fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]];
222               rtmp=factor2 * (m == 0 ? 1 : m);
223               fctp->constraint_weights[i]=rtmp;
224               fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]]=m+rtmp;
225            } else {
226               fctp->constraint_weights[i]=0;
227            }
228       }
229       fctp->constraint_factor=factor;
230    }
231    
232    void Paso_FCTransportProblem_insertConstraint(Paso_FCTransportProblem* fctp,  const double* r,  double* source)
233    {
234       dim_t i, n;
235       n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
236    
237       #pragma omp for schedule(static) private(i,m,rtmp)
238       for (i=0;i<n;++i) source[i]+=fctp->constraint_weights[i] * r[i];
239    }

Legend:
Removed from v.2196  
changed lines
  Added in v.2197

  ViewVC Help
Powered by ViewVC 1.1.26