Diff of /trunk/paso/src/FCT_Solver.cpp

revision 4846 by caltinay, Tue Apr 8 23:55:38 2014 UTC revision 4852 by caltinay, Wed Apr 9 03:58:08 2014 UTC
# Line 233  err_t Paso_FCT_Solver_update_LCN(Paso_FC Line 233  err_t Paso_FCT_Solver_update_LCN(Paso_FC
235
236     Paso_Scale(n, b,fct_solver->omega );     paso::util::scale(n, b,fct_solver->omega );
237     /* solve (m-dt/2*L) u = b in the form (omega*m-L) u = b * omega with omega*dt/2=1 */     /* solve (m-dt/2*L) u = b in the form (omega*m-L) u = b * omega with omega*dt/2=1 */
238     #pragma omp for private(i) schedule(static)     #pragma omp for private(i) schedule(static)
239     for (i = 0; i < n; ++i) {     for (i = 0; i < n; ++i) {
# Line 243  err_t Paso_FCT_Solver_update_LCN(Paso_FC Line 243  err_t Paso_FCT_Solver_update_LCN(Paso_FC
243         }         }
244     }     }
245     /* initial guess is u<- -u + 2*u_tilde */     /* initial guess is u<- -u + 2*u_tilde */
246     Paso_Update(n, -1., u, 2., fct_solver->flux_limiter->u_tilde);     paso::util::update(n, -1., u, 2., fct_solver->flux_limiter->u_tilde);
247
248     sweep_max = MAX((int) (- 2 * log(RTOL)/log(2.)-0.5),1);     sweep_max = MAX((int) (- 2 * log(RTOL)/log(2.)-0.5),1);
249     norm_u_tilde=Paso_lsup(n, fct_solver->flux_limiter->u_tilde, fct_solver->flux_limiter->mpi_info);     norm_u_tilde=paso::util::lsup(n, fct_solver->flux_limiter->u_tilde, fct_solver->flux_limiter->mpi_info);
250     if (options->verbose) {     if (options->verbose) {
251         printf("Paso_FCT_Solver_update_LCN: u_tilda lsup = %e (rtol = %e, max. sweeps = %d)\n",norm_u_tilde,RTOL * norm_u_tilde ,sweep_max);         printf("Paso_FCT_Solver_update_LCN: u_tilda lsup = %e (rtol = %e, max. sweeps = %d)\n",norm_u_tilde,RTOL * norm_u_tilde ,sweep_max);
252     }     }
# Line 311  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 311  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
311      /* u_tilda_connector is completed */      /* u_tilda_connector is completed */
312      /************************************************************************/      /************************************************************************/
313      /* calculate stopping criterion */      /* calculate stopping criterion */
314      norm_u_tilde=Paso_lsup(n, flux_limiter->u_tilde, flux_limiter->mpi_info);      norm_u_tilde=paso::util::lsup(n, flux_limiter->u_tilde, flux_limiter->mpi_info);
315      ATOL= rtol * norm_u_tilde + atol ;      ATOL= rtol * norm_u_tilde + atol ;
316      if (options->verbose) printf("Paso_FCT_Solver_updateNL: iteration starts u_tilda lsup = %e (abs. tol = %e)\n",norm_u_tilde,ATOL);      if (options->verbose) printf("Paso_FCT_Solver_updateNL: iteration starts u_tilda lsup = %e (abs. tol = %e)\n",norm_u_tilde,ATOL);
317      /* ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// */      /* ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// */
318
319      /* u_old is an initial guess for u*/      /* u_old is an initial guess for u*/
320      Paso_Copy(n,u,u_old);      paso::util::copy(n,u,u_old);
321
322      while ( (!converged) && (!diverged) && (! max_m_reached) && Esys_noError()) {      while ( (!converged) && (!diverged) && (! max_m_reached) && Esys_noError()) {
323          fct_solver->u_coupler->startCollect(u);          fct_solver->u_coupler->startCollect(u);
# Line 347  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 347  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
347                Paso_FCT_Solver_setMuPaLu(z, fctp->lumped_mass_matrix, fct_solver->u_coupler, dt/2, fctp->iteration_matrix);                Paso_FCT_Solver_setMuPaLu(z, fctp->lumped_mass_matrix, fct_solver->u_coupler, dt/2, fctp->iteration_matrix);
348            }            }
349
350            Paso_Update(n,-1.,z,1.,b);  /* z=b-z */            paso::util::update(n,-1.,z,1.,b);  /* z=b-z */
351
352            /* z_i += sum_{j} limitation factor_{ij} * antidiffusive_flux_{ij} */            /* z_i += sum_{j} limitation factor_{ij} * antidiffusive_flux_{ij} */
# Line 355  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 355  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
355            /* we solve (m/omega - L ) * du = z */            /* we solve (m/omega - L ) * du = z */
356            if (fct_solver->method == PASO_BACKWARD_EULER) {            if (fct_solver->method == PASO_BACKWARD_EULER) {
357                dim_t cntIter = options->iter_max;                dim_t cntIter = options->iter_max;
358                double tol= Paso_l2(n, z, fctp->mpi_info) ;                double tol= paso::util::l2(n, z, fctp->mpi_info) ;
359
360                if ( m ==0) {                if ( m ==0) {
361                    tol *=0.5;                    tol *=0.5;
# Line 363  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 363  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
363                    tol *= MIN(MAX(rate*rate, 1e-2), 0.5);                    tol *= MIN(MAX(rate*rate, 1e-2), 0.5);
364                }                }
365                /* use BiCGSTab with jacobi preconditioner ( m - omega * L ) */                /* use BiCGSTab with jacobi preconditioner ( m - omega * L ) */
366                Paso_zeroes(n,du);                paso::util::zeroes(n,du);
367                errorCode = Paso_Solver_BiCGStab(fctp->iteration_matrix, z, du, &cntIter, &tol, pp);                errorCode = Paso_Solver_BiCGStab(fctp->iteration_matrix, z, du, &cntIter, &tol, pp);
368
369                /* errorCode =  Paso_Solver_GMRES(fctp->iteration_matrix, z, du, &cntIter, &tol, 10, 2000, pp); */                /* errorCode =  Paso_Solver_GMRES(fctp->iteration_matrix, z, du, &cntIter, &tol, 10, 2000, pp); */
# Line 379  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 379  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
379                options->num_iter++;                options->num_iter++;
380             }             }
381
382             Paso_Update(n,1.,u,fct_solver->omega,du);            paso::util::update(n,1.,u,fct_solver->omega,du);
383             norm_du_old=norm_du;             norm_du_old=norm_du;
384             norm_du=Paso_lsup(n,du, fctp->mpi_info);             norm_du=paso::util::lsup(n,du, fctp->mpi_info);
385             if (m ==0) {             if (m ==0) {
386                  if (options->verbose) printf("Paso_FCT_Solver_updateNL: step %d: increment= %e\n",m+1, norm_du * fct_solver->omega);                  if (options->verbose) printf("Paso_FCT_Solver_updateNL: step %d: increment= %e\n",m+1, norm_du * fct_solver->omega);
387             } else {             } else {

Legend:
 Removed from v.4846 changed lines Added in v.4852