/[escript]/trunk/paso/src/SolverFCT_solve.c
ViewVC logotype

Diff of /trunk/paso/src/SolverFCT_solve.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2196 by phornby, Thu Nov 20 16:10:10 2008 UTC revision 2197 by gross, Thu Jan 8 05:49:16 2009 UTC
# Line 181  void Paso_SolverFCT_solve(Paso_FCTranspo Line 181  void Paso_SolverFCT_solve(Paso_FCTranspo
181          } else {          } else {
182               dt2=dt;               dt2=dt;
183          }          }
184          while(t<dt && Paso_noError()) {          while( (t<dt*(1.-sqrt(EPSILON))) && Paso_noError()) {
185              printf("substep step %d at t=%e (step size= %e)\n",i_substeps+1,t+dt2,dt2);              printf("substep step %d at t=%e (step size= %e)\n",i_substeps+1,t+dt2,dt2);
186              Paso_FCT_setUp(fctp,dt2,sourceN,sourceP,b_n,uTilde_n,uTilde_n_coupler,QN_n,QN_n_coupler,QP_n,QP_n_coupler,              Paso_FCT_setUp(fctp,dt2,sourceN,sourceP,b_n,uTilde_n,uTilde_n_coupler,QN_n,QN_n_coupler,QP_n,QP_n_coupler,
187                             options,&pp);                             options,&pp);
# Line 198  void Paso_SolverFCT_solve(Paso_FCTranspo Line 198  void Paso_SolverFCT_solve(Paso_FCTranspo
198                       * now we solve the linear system to get the correction dt:                       * now we solve the linear system to get the correction dt:
199                       *                       *
200                       */                       */
201    /* for (i=0;i<n;++i) printf("%d %f \n", i,z_m[i]); */
202                       if (fctp->theta > 0) {                       if (fctp->theta > 0) {
203                            omega=1./(dt2*fctp->theta);                            omega=1./(dt2*fctp->theta);
204                            Paso_Solver_solvePreconditioner(fctp->iteration_matrix,du_m,z_m);                              Paso_Solver_solvePreconditioner(fctp->iteration_matrix,du_m,z_m);  
# Line 278  double Paso_FCTransportProblem_getSafeTi Line 279  double Paso_FCTransportProblem_getSafeTi
279     double dt_max, dt_max_loc;     double dt_max, dt_max_loc;
280     register double l_ii,m;     register double l_ii,m;
281     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
282     if (! fctp->valid_matrices) {     if ( ! fctp->valid_matrices) {
283          fctp->dt_max=LARGE_POSITIVE_FLOAT;        /* extract the row sum of the advective part */
284          /* extract the row sum of the advective part */        Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);
         Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);  
285    
286          /* set low order transport operator */        /* set low order transport operator */
287          Paso_FCTransportProblem_setLowOrderOperator(fctp);        Paso_FCTransportProblem_setLowOrderOperator(fctp);
288    
289          if (Paso_noError()) {
290          /*          /*
291           *  calculate time step size:                                                     *  calculate time step size:                                          
292          */          */
293          dt_max=LARGE_POSITIVE_FLOAT;          dt_max=LARGE_POSITIVE_FLOAT;
294          if (fctp->theta < 1.) {          #pragma omp parallel private(dt_max_loc)
295              #pragma omp parallel private(dt_max_loc)          {
             {  
296                 dt_max_loc=LARGE_POSITIVE_FLOAT;                 dt_max_loc=LARGE_POSITIVE_FLOAT;
297                 #pragma omp for schedule(static) private(i,l_ii,m)                 #pragma omp for schedule(static) private(i,l_ii,m)
298                 for (i=0;i<n;++i) {                 for (i=0;i<n;++i) {
# Line 305  double Paso_FCTransportProblem_getSafeTi Line 306  double Paso_FCTransportProblem_getSafeTi
306                 {                 {
307                    dt_max=MIN(dt_max,dt_max_loc);                    dt_max=MIN(dt_max,dt_max_loc);
308                 }                 }
             }  
             #ifdef PASO_MPI  
                dt_max_loc = dt_max;  
                MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);  
             #endif  
             if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=1./(1.-fctp->theta);  
309           }           }
310             #ifdef PASO_MPI
311                 dt_max_loc = dt_max;
312                 MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
313             #endif
314    printf("dt_max = %e %e\n",dt_max, fctp->dt_factor);
315             if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=fctp->dt_factor;
316           if (dt_max <= 0.)  {           if (dt_max <= 0.)  {
317              Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");              Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
318           } else {           } else {
319              if (dt_max<LARGE_POSITIVE_FLOAT)              if (dt_max<LARGE_POSITIVE_FLOAT)
320                 printf("maximum time step size is %e (theta = %e).\n",dt_max,fctp->theta);                   printf("maximum time step size is %e (theta = %e).\n",dt_max,fctp->theta);  
321          }           }
322          fctp->dt_max=dt_max;           fctp->dt_max=dt_max;
323          fctp->valid_matrices=Paso_noError();           fctp->valid_matrices=Paso_noError();
324          }
325     }     }
326     return fctp->dt_max;     return fctp->dt_max;
327  }  }
Legend:
Removed from v.2196  
changed lines
  Added in v.2197
  ViewVC Help
Powered by ViewVC 1.1.26