/[escript]/trunk/paso/src/FCT_Solver.cpp
ViewVC logotype

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

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

revision 4810 by caltinay, Thu Mar 27 07:00:30 2014 UTC revision 4817 by caltinay, Fri Mar 28 08:04:09 2014 UTC
# Line 53  Paso_FCT_Solver* Paso_FCT_Solver_alloc(P Line 53  Paso_FCT_Solver* Paso_FCT_Solver_alloc(P
53           out->du = NULL;           out->du = NULL;
54           out->z=NULL;           out->z=NULL;
55          }          }
56      out->u_coupler = paso::Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp), blockSize);      out->u_coupler.reset(new paso::Coupler(Paso_TransportProblem_borrowConnector(fctp), blockSize));
57          out->u_old_coupler = paso::Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp), blockSize);      out->u_old_coupler.reset(new paso::Coupler(Paso_TransportProblem_borrowConnector(fctp), blockSize));
58      out->omega=0;      out->omega=0;
59            
60      if ( options->ode_solver == PASO_LINEAR_CRANK_NICOLSON ) {      if ( options->ode_solver == PASO_LINEAR_CRANK_NICOLSON ) {
61           out->method   = PASO_LINEAR_CRANK_NICOLSON;           out->method   = PASO_LINEAR_CRANK_NICOLSON;
# Line 80  Paso_FCT_Solver* Paso_FCT_Solver_alloc(P Line 80  Paso_FCT_Solver* Paso_FCT_Solver_alloc(P
80  void Paso_FCT_Solver_free(Paso_FCT_Solver *in)  void Paso_FCT_Solver_free(Paso_FCT_Solver *in)
81  {  {
82      if (in != NULL) {      if (in != NULL) {
83            Paso_TransportProblem_free(in->transportproblem);          Paso_TransportProblem_free(in->transportproblem);
84        Paso_FCT_FluxLimiter_free(in->flux_limiter);          Paso_FCT_FluxLimiter_free(in->flux_limiter);
85            Esys_MPIInfo_free(in->mpi_info);          Esys_MPIInfo_free(in->mpi_info);
86            paso::Coupler_free(in->u_old_coupler);          delete[] in->b;
87            paso::Coupler_free(in->u_coupler);          delete[] in->z;
88            delete[] in->du;
89        delete[] in->b;          delete in;
       delete[] in->z;  
       delete[] in->du;  
           delete in;  
         
90      }      }
91  }  }
92    
# Line 218  err_t Paso_FCT_Solver_update_LCN(Paso_FC Line 214  err_t Paso_FCT_Solver_update_LCN(Paso_FC
214      err_t errorCode = SOLVER_NO_ERROR;      err_t errorCode = SOLVER_NO_ERROR;
215      double norm_u_tilde;      double norm_u_tilde;
216        
217      paso::Coupler_startCollect(fct_solver->u_old_coupler,u_old);      fct_solver->u_old_coupler->startCollect(u_old);
218      paso::Coupler_finishCollect(fct_solver->u_old_coupler);      fct_solver->u_old_coupler->finishCollect();
219            
220      /* b[i]=m*u_tilde[i] = m u_old[i] + dt/2 sum_{j <> i} l_{ij}*(u_old[j]-u_old[i])        /* b[i]=m*u_tilde[i] = m u_old[i] + dt/2 sum_{j <> i} l_{ij}*(u_old[j]-u_old[i])  
221             = u_tilde[i]   = u_old[i] where constraint m<0.             = u_tilde[i]   = u_old[i] where constraint m<0.
# Line 297  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 293  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
293     bool converged=FALSE, max_m_reached=FALSE,diverged=FALSE;       bool converged=FALSE, max_m_reached=FALSE,diverged=FALSE;  
294     options->num_iter=0;     options->num_iter=0;
295        
296     paso::Coupler_startCollect(fct_solver->u_old_coupler,u_old);     fct_solver->u_old_coupler->startCollect(u_old);
297     paso::Coupler_finishCollect(fct_solver->u_old_coupler);     fct_solver->u_old_coupler->finishCollect();
298      /* prepare u_tilda and flux limiter */      /* prepare u_tilda and flux limiter */
299      if ( fct_solver->method == PASO_BACKWARD_EULER ) {      if ( fct_solver->method == PASO_BACKWARD_EULER ) {
300            /* b[i]=m_i* u_old[i] */            /* b[i]=m_i* u_old[i] */
# Line 330  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 326  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
326      Paso_Copy(n,u,u_old);      Paso_Copy(n,u,u_old);
327            
328      while ( (!converged) && (!diverged) && (! max_m_reached) && Esys_noError()) {      while ( (!converged) && (!diverged) && (! max_m_reached) && Esys_noError()) {
329          paso::Coupler_startCollect(fct_solver->u_coupler,u);          fct_solver->u_coupler->startCollect(u);
330          paso::Coupler_finishCollect(fct_solver->u_coupler);          fct_solver->u_coupler->finishCollect();
331    
332           /*  set antidiffusive_flux_{ij} for u */           /*  set antidiffusive_flux_{ij} for u */
333           if (fct_solver->method == PASO_BACKWARD_EULER) {           if (fct_solver->method == PASO_BACKWARD_EULER) {
# Line 443  err_t Paso_FCT_Solver_updateNL(Paso_FCT_ Line 439  err_t Paso_FCT_Solver_updateNL(Paso_FCT_
439    
440  void Paso_FCT_setAntiDiffusionFlux_CN(Paso_SystemMatrix *flux_matrix,  void Paso_FCT_setAntiDiffusionFlux_CN(Paso_SystemMatrix *flux_matrix,
441                                        const Paso_TransportProblem* fct,                                        const Paso_TransportProblem* fct,
442                        const double dt,                                        const double dt,
443                            const paso::Coupler* u_coupler,                                          paso::const_Coupler_ptr u_coupler,  
444                        const paso::Coupler* u_old_coupler)                                        paso::const_Coupler_ptr u_old_coupler)
445  {  {
446    dim_t i;    dim_t i;
447    index_t iptr_ij;    index_t iptr_ij;
448        
449    const double *u    = paso::Coupler_borrowLocalData(u_coupler);    const double *u    = u_coupler->borrowLocalData();
450    const double *u_old= paso::Coupler_borrowLocalData(u_old_coupler);    const double *u_old= u_old_coupler->borrowLocalData();
451    const double *remote_u=paso::Coupler_borrowRemoteData(u_coupler);    const double *remote_u=u_coupler->borrowRemoteData();
452    const double *remote_u_old=paso::Coupler_borrowRemoteData(u_old_coupler);    const double *remote_u_old=u_old_coupler->borrowRemoteData();
453    const double dt_half= dt/2;    const double dt_half= dt/2;
454    const paso::SystemMatrixPattern *pattern=fct->iteration_matrix->pattern;    const paso::SystemMatrixPattern* pattern=fct->iteration_matrix->pattern;
455    const dim_t n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);    const dim_t n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);
456    
457    #pragma omp parallel for schedule(static) private(i, iptr_ij)    #pragma omp parallel for schedule(static) private(i, iptr_ij)
# Line 489  void Paso_FCT_setAntiDiffusionFlux_CN(Pa Line 485  void Paso_FCT_setAntiDiffusionFlux_CN(Pa
485    
486  void Paso_FCT_setAntiDiffusionFlux_BE(Paso_SystemMatrix *flux_matrix,  void Paso_FCT_setAntiDiffusionFlux_BE(Paso_SystemMatrix *flux_matrix,
487                                        const Paso_TransportProblem* fct,                                        const Paso_TransportProblem* fct,
488                        const double dt,                                        const double dt,
489                            const paso::Coupler* u_coupler,                                          paso::const_Coupler_ptr u_coupler,  
490                        const paso::Coupler* u_old_coupler)                                        paso::const_Coupler_ptr u_old_coupler)
491  {  {
492    dim_t i;    dim_t i;
493    index_t iptr_ij;    index_t iptr_ij;
494        
495    const double *u=paso::Coupler_borrowLocalData(u_coupler);    const double *u = u_coupler->borrowLocalData();
496    const double *u_old= paso::Coupler_borrowLocalData(u_old_coupler);    const double *u_old = u_old_coupler->borrowLocalData();
497    const double *remote_u=paso::Coupler_borrowRemoteData(u_coupler);    const double *remote_u = u_coupler->borrowRemoteData();
498    const double *remote_u_old=paso::Coupler_borrowRemoteData(u_old_coupler);    const double *remote_u_old = u_old_coupler->borrowRemoteData();
499    const paso::SystemMatrixPattern *pattern=fct->iteration_matrix->pattern;    const paso::SystemMatrixPattern *pattern=fct->iteration_matrix->pattern;
500    const dim_t  n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);    const dim_t  n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);
501    
# Line 542  void Paso_FCT_setAntiDiffusionFlux_BE(Pa Line 538  void Paso_FCT_setAntiDiffusionFlux_BE(Pa
538   */   */
539        
540  void Paso_FCT_setAntiDiffusionFlux_linearCN(Paso_SystemMatrix *flux_matrix,  void Paso_FCT_setAntiDiffusionFlux_linearCN(Paso_SystemMatrix *flux_matrix,
541                                              const Paso_TransportProblem* fct,                          const Paso_TransportProblem* fct, const double dt,
542                              const double dt,                          paso::const_Coupler_ptr u_tilde_coupler,
543                                  const paso::Coupler* u_tilde_coupler,                            paso::const_Coupler_ptr u_old_coupler)
                             const paso::Coupler* u_old_coupler)  
544  {  {
545    dim_t i;    dim_t i;
546    index_t iptr_ij;    index_t iptr_ij;
547        
548    const double *u_tilde=paso::Coupler_borrowLocalData(u_tilde_coupler);    const double *u_tilde = u_tilde_coupler->borrowLocalData();
549    const double *u_old= paso::Coupler_borrowLocalData(u_old_coupler);    const double *u_old = u_old_coupler->borrowLocalData();
550    const double *remote_u_tilde=paso::Coupler_borrowRemoteData(u_tilde_coupler);    const double *remote_u_tilde = u_tilde_coupler->borrowRemoteData();
551    const double *remote_u_old=paso::Coupler_borrowRemoteData(u_old_coupler);    const double *remote_u_old = u_old_coupler->borrowRemoteData();
552    const paso::SystemMatrixPattern *pattern=fct->iteration_matrix->pattern;    const paso::SystemMatrixPattern *pattern=fct->iteration_matrix->pattern;
553    const dim_t n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);    const dim_t n=Paso_SystemMatrix_getTotalNumRows(fct->iteration_matrix);
554    
# Line 684  printf("a[%d,%d]=%e\n",j,i,rtmp2); Line 679  printf("a[%d,%d]=%e\n",j,i,rtmp2);
679   *   *
680   */   */
681  void Paso_FCT_Solver_setMuPaLu(double* out,  void Paso_FCT_Solver_setMuPaLu(double* out,
682                                const double* M,                                 const double* M,
683                                const paso::Coupler* u_coupler,                                 paso::const_Coupler_ptr u_coupler,
684                                const double a,                                 const double a,
685                                const Paso_SystemMatrix *L)                                 const Paso_SystemMatrix *L)
686  {  {
687    dim_t i;    dim_t i;
688    const paso::SystemMatrixPattern *pattern = L->pattern;    const paso::SystemMatrixPattern *pattern = L->pattern;
689    const double *u=paso::Coupler_borrowLocalData(u_coupler);    const double *u = u_coupler->borrowLocalData();
690    const double *remote_u=paso::Coupler_borrowRemoteData(u_coupler);    const double *remote_u = u_coupler->borrowRemoteData();
691    index_t iptr_ij;    index_t iptr_ij;
692    const dim_t n=Paso_SystemMatrix_getTotalNumRows(L);    const dim_t n=Paso_SystemMatrix_getTotalNumRows(L);
693    

Legend:
Removed from v.4810  
changed lines
  Added in v.4817

  ViewVC Help
Powered by ViewVC 1.1.26