/[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 1916 by ksteube, Fri Sep 26 00:19:18 2008 UTC revision 1917 by phornby, Thu Oct 23 09:09:31 2008 UTC
# Line 69  err_t Paso_FCT_setUpRightHandSide(Paso_F Line 69  err_t Paso_FCT_setUpRightHandSide(Paso_F
69                                    double *RN_m, Paso_Coupler* RN_m_coupler, double* RP_m, Paso_Coupler* RP_m_coupler, const double *sourceN,                                    double *RN_m, Paso_Coupler* RN_m_coupler, double* RP_m, Paso_Coupler* RP_m_coupler, const double *sourceN,
70                                    Paso_Performance* pp)                                    Paso_Performance* pp)
71  {  {
    dim_t i;  
72     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
    register double m, rtmp;  
73     /* distribute u */     /* distribute u */
74     Paso_Coupler_startCollect(u_m_coupler,u_m);     Paso_Coupler_startCollect(u_m_coupler,u_m);
75     Paso_Coupler_finishCollect(u_m_coupler);     Paso_Coupler_finishCollect(u_m_coupler);
# Line 111  err_t Paso_FCT_setUpRightHandSide(Paso_F Line 109  err_t Paso_FCT_setUpRightHandSide(Paso_F
109    
110  void Paso_SolverFCT_solve(Paso_FCTransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options) {  void Paso_SolverFCT_solve(Paso_FCTransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options) {
111     const dim_t FAILURES_MAX=5;     const dim_t FAILURES_MAX=5;
112     err_t error_code;     dim_t m, i_substeps, Failed, i;
    dim_t m,n_substeps, i_substeps, Failed, i, iter;  
    err_t errorCode;  
113     double *z_m=NULL, *b_n=NULL, *sourceP=NULL, *sourceN=NULL, *uTilde_n=NULL, *QN_n=NULL, *QP_n=NULL, *RN_m=NULL, *RP_m=NULL, *du_m=NULL;     double *z_m=NULL, *b_n=NULL, *sourceP=NULL, *sourceN=NULL, *uTilde_n=NULL, *QN_n=NULL, *QP_n=NULL, *RN_m=NULL, *RP_m=NULL, *du_m=NULL;
114     Paso_Coupler *QN_n_coupler=NULL, *QP_n_coupler=NULL, *RN_m_coupler=NULL, *RP_m_coupler=NULL, *uTilde_n_coupler=NULL, *u_m_coupler=NULL;     Paso_Coupler *QN_n_coupler=NULL, *QP_n_coupler=NULL, *RN_m_coupler=NULL, *RP_m_coupler=NULL, *uTilde_n_coupler=NULL, *u_m_coupler=NULL;
115     Paso_SystemMatrix *flux_matrix_m=NULL;     Paso_SystemMatrix *flux_matrix_m=NULL;
116     double dt_max, dt2,t, norm_u_m, omega, norm_du_m, tol;     double dt_max, dt2,t, norm_u_m, omega, norm_du_m;
117     register double mass, rtmp;     register double mass, rtmp;
118     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
119     dim_t blockSize=Paso_FCTransportProblem_getBlockSize(fctp);     dim_t blockSize=Paso_FCTransportProblem_getBlockSize(fctp);
# Line 299  double Paso_FCTransportProblem_getSafeTi Line 295  double Paso_FCTransportProblem_getSafeTi
295                 dt_max_loc=LARGE_POSITIVE_FLOAT;                 dt_max_loc=LARGE_POSITIVE_FLOAT;
296                 #pragma omp for schedule(static) private(i,l_ii,m)                 #pragma omp for schedule(static) private(i,l_ii,m)
297                 for (i=0;i<n;++i) {                 for (i=0;i<n;++i) {
298                      l_ii=fctp->main_diagonal_low_order_transport_matrix[i];                    l_ii=fctp->main_diagonal_low_order_transport_matrix[i];
299                      m=fctp->lumped_mass_matrix[i];                    m=fctp->lumped_mass_matrix[i];
300                      if ( (l_ii<0 && m>0.) || (l_ii>0 && m<0) ) {                    if ( (l_ii<0 && m>0.) || (l_ii>0 && m<0) ) {
301                          dt_max_loc=MIN(dt_max_loc,-m/l_ii);                       dt_max_loc=MIN(dt_max_loc,-m/l_ii);
302                      }                    }
303                  }                 }
304                  #pragma omp critical                  #pragma omp critical
305                  {                 {
306                      dt_max=MIN(dt_max,dt_max_loc);                    dt_max=MIN(dt_max,dt_max_loc);
307                  }                 }
308              }              }
309              #ifdef PASO_MPI              #ifdef PASO_MPI
310                 dt_max_loc = dt_max;                 dt_max_loc = dt_max;
311             MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);                 MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
312               #endif              #endif
313               if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=1./(1.-fctp->theta);              if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=1./(1.-fctp->theta);
314           }           }
315           if (dt_max <= 0.)  {           if (dt_max <= 0.)  {
316              Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");              Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
317           } else {           } else {
318             if (dt_max<LARGE_POSITIVE_FLOAT) printf("maximum time step size is %e (theta = %e).\n",dt_max,fctp->theta);                if (dt_max<LARGE_POSITIVE_FLOAT)
319                   printf("maximum time step size is %e (theta = %e).\n",dt_max,fctp->theta);  
320          }          }
321          fctp->dt_max=dt_max;          fctp->dt_max=dt_max;
322          fctp->valid_matrices=Paso_noError();          fctp->valid_matrices=Paso_noError();
Legend:
Removed from v.1916  
changed lines
  Added in v.1917
  ViewVC Help
Powered by ViewVC 1.1.26