/[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 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC revision 1812 by ksteube, Fri Sep 26 00:19:18 2008 UTC
# Line 40  Line 40 
40  #include <mpi.h>  #include <mpi.h>
41  #endif  #endif
42    
43    /*
44     * inserts the source term into the problem
45     */
46    void Paso_FCT_setSource(const dim_t n,const double *source, double* sourceN, double* sourceP)
47    {
48       dim_t i;
49       register double rtmp;
50       /*
51        * seperate source into positive and negative part:
52        */
53       #pragma omp parallel for private(i,rtmp)
54       for (i = 0; i < n; ++i) {
55           rtmp=source[i];
56           if (rtmp <0) {
57              sourceN[i]=-rtmp;
58              sourceP[i]=0;
59           } else {
60              sourceN[i]= 0;
61              sourceP[i]= rtmp;
62           }
63       }
64    }
65    
66    err_t Paso_FCT_setUpRightHandSide(Paso_FCTransportProblem* fctp, const double dt, const double *u_m, Paso_Coupler* u_m_coupler,  double * z_m,
67                                      Paso_SystemMatrix* flux_matrix, Paso_Coupler* uTilde_coupler, const double *b,
68                                      Paso_Coupler* QN_coupler, Paso_Coupler* QP_coupler,
69                                      double *RN_m, Paso_Coupler* RN_m_coupler, double* RP_m, Paso_Coupler* RP_m_coupler, const double *sourceN,
70                                      Paso_Performance* pp)
71    {
72       dim_t i;
73       const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
74       register double m, rtmp;
75       /* distribute u */
76       Paso_Coupler_startCollect(u_m_coupler,u_m);
77       Paso_Coupler_finishCollect(u_m_coupler);
78       /*
79        *  set the ant diffusion fluxes:
80        *
81        */
82       Paso_FCTransportProblem_setAntiDiffusionFlux(dt,fctp,flux_matrix,u_m_coupler);
83       /*
84        *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(\tilde{u}[i]- \tilde{u}[j])<=0
85        *
86        */
87       Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(flux_matrix,uTilde_coupler);
88       /*
89        *  set flux limiters RN_m,RP_m
90        *
91        */
92       Paso_FCTransportProblem_setRs(flux_matrix,fctp->lumped_mass_matrix,QN_coupler,QP_coupler,RN_m,RP_m);
93       Paso_Coupler_startCollect(RN_m_coupler,RN_m);
94       Paso_Coupler_startCollect(RP_m_coupler,RP_m);
95        /*
96         * z_m[i]=b[i] - (m_i*u_m[i] - dt*theta*sum_{j<>i} l_{ij} (u_m[j]-u_m[i]) + dt q^-[i])
97         *
98         * note that iteration_matrix stores the negative values of the
99         * low order transport matrix l therefore a=dt*fctp->theta is used.
100         */
101       Paso_SolverFCT_setMuPaLuPbQ(z_m,fctp->lumped_mass_matrix, u_m_coupler,dt*fctp->theta,fctp->iteration_matrix,dt,sourceN);
102       /* z_m=b-z_m */
103       Paso_Update(n,-1.,z_m,1.,b);
104    
105       Paso_Coupler_finishCollect(RN_m_coupler);
106       Paso_Coupler_finishCollect(RP_m_coupler);
107       /* add corrected fluxes into z_m */
108       Paso_FCTransportProblem_addCorrectedFluxes(z_m,flux_matrix,RN_m_coupler,RP_m_coupler);
109       return NO_ERROR;
110    }
111    
112  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) {
113     const dim_t FAILURES_MAX=5;     const dim_t FAILURES_MAX=5;
114     err_t error_code;     err_t error_code;
# Line 117  void Paso_SolverFCT_solve(Paso_FCTranspo Line 186  void Paso_SolverFCT_solve(Paso_FCTranspo
186               dt2=dt;               dt2=dt;
187          }          }
188          while(t<dt && Paso_noError()) {          while(t<dt && Paso_noError()) {
189              printf("substep step %ld 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);
190              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,
191                             options,&pp);                             options,&pp);
192              /* now the iteration starts */              /* now the iteration starts */
# Line 153  void Paso_SolverFCT_solve(Paso_FCTranspo Line 222  void Paso_SolverFCT_solve(Paso_FCTranspo
222                     }                     }
223                     norm_u_m=Paso_lsup(n,u, fctp->mpi_info);                     norm_u_m=Paso_lsup(n,u, fctp->mpi_info);
224                     norm_du_m=Paso_lsup(n,du_m, fctp->mpi_info)*omega;                     norm_du_m=Paso_lsup(n,du_m, fctp->mpi_info)*omega;
225                     printf("iteration step %ld completed: norm increment= %e (tolerance = %e)\n",m+1, norm_du_m, rtol * norm_u_m + atol);                     printf("iteration step %d completed: norm increment= %e (tolerance = %e)\n",m+1, norm_du_m, rtol * norm_u_m + atol);
226    
227                     max_m_reached=(m>max_m);                     max_m_reached=(m>max_m);
228                     converged=(norm_du_m <= rtol * norm_u_m + atol);                     converged=(norm_du_m <= rtol * norm_u_m + atol);
# Line 161  void Paso_SolverFCT_solve(Paso_FCTranspo Line 230  void Paso_SolverFCT_solve(Paso_FCTranspo
230              }              }
231              if (converged) {              if (converged) {
232                      Failed=0;                      Failed=0;
233                      #pragma omp parallel for schedule(static) private(i)                      /* #pragma omp parallel for schedule(static) private(i) */
234                      Paso_Copy(n,fctp->u,u);                      Paso_Copy(n,fctp->u,u);
235                      i_substeps++;                      i_substeps++;
236                      t+=dt2;                      t+=dt2;
# Line 257  double Paso_FCTransportProblem_getSafeTi Line 326  double Paso_FCTransportProblem_getSafeTi
326     }     }
327     return fctp->dt_max;     return fctp->dt_max;
328  }  }
329  /*  void Paso_FCT_setUp(Paso_FCTransportProblem* fctp, const double dt, const double *sourceN, const double *sourceP, double* b, double* uTilde,
  * inserts the source term into the problem  
  */  
 err_t Paso_FCT_setSource(const dim_t n,const double *source, double* sourceN, double* sourceP)  
 {  
    dim_t i;  
    register double rtmp;  
    /*  
     * seperate source into positive and negative part:  
     */  
    #pragma omp parallel for private(i,rtmp)  
    for (i = 0; i < n; ++i) {  
        rtmp=source[i];  
        if (rtmp <0) {  
           sourceN[i]=-rtmp;  
           sourceP[i]=0;  
        } else {  
           sourceN[i]= 0;  
           sourceP[i]= rtmp;  
        }  
    }  
 }  
 err_t Paso_FCT_setUp(Paso_FCTransportProblem* fctp, const double dt, const double *sourceN, const double *sourceP, double* b, double* uTilde,  
330                       Paso_Coupler* uTilde_coupler, double *QN, Paso_Coupler* QN_coupler, double *QP, Paso_Coupler* QP_coupler,                       Paso_Coupler* uTilde_coupler, double *QN, Paso_Coupler* QN_coupler, double *QP, Paso_Coupler* QP_coupler,
331                       Paso_Options* options, Paso_Performance* pp)                       Paso_Options* options, Paso_Performance* pp)
332  {  {
# Line 353  err_t Paso_FCT_setUp(Paso_FCTransportPro Line 400  err_t Paso_FCT_setUp(Paso_FCTransportPro
400       Paso_Coupler_finishCollect(QN_coupler);       Paso_Coupler_finishCollect(QN_coupler);
401       Paso_Coupler_finishCollect(QP_coupler);       Paso_Coupler_finishCollect(QP_coupler);
402  }  }
 err_t Paso_FCT_setUpRightHandSide(Paso_FCTransportProblem* fctp, const double dt, const double *u_m, Paso_Coupler* u_m_coupler,  double * z_m,  
                                   Paso_SystemMatrix* flux_matrix, Paso_Coupler* uTilde_coupler, const double *b,  
                                   Paso_Coupler* QN_coupler, Paso_Coupler* QP_coupler,  
                                   double *RN_m, Paso_Coupler* RN_m_coupler, double* RP_m, Paso_Coupler* RP_m_coupler, const double *sourceN,  
                                   Paso_Performance* pp)  
 {  
    dim_t i;  
    const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);  
    register double m, rtmp;  
    /* distribute u */  
    Paso_Coupler_startCollect(u_m_coupler,u_m);  
    Paso_Coupler_finishCollect(u_m_coupler);  
    /*  
     *  set the ant diffusion fluxes:  
     *  
     */  
    Paso_FCTransportProblem_setAntiDiffusionFlux(dt,fctp,flux_matrix,u_m_coupler);  
    /*  
     *  apply pre flux-correction: f_{ij}:=0 if f_{ij}*(\tilde{u}[i]- \tilde{u}[j])<=0  
     *  
     */  
    Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(flux_matrix,uTilde_coupler);  
    /*  
     *  set flux limiters RN_m,RP_m  
     *  
     */  
    Paso_FCTransportProblem_setRs(flux_matrix,fctp->lumped_mass_matrix,QN_coupler,QP_coupler,RN_m,RP_m);  
    Paso_Coupler_startCollect(RN_m_coupler,RN_m);  
    Paso_Coupler_startCollect(RP_m_coupler,RP_m);  
     /*  
      * z_m[i]=b[i] - (m_i*u_m[i] - dt*theta*sum_{j<>i} l_{ij} (u_m[j]-u_m[i]) + dt q^-[i])  
      *  
      * note that iteration_matrix stores the negative values of the  
      * low order transport matrix l therefore a=dt*fctp->theta is used.  
      */  
    Paso_SolverFCT_setMuPaLuPbQ(z_m,fctp->lumped_mass_matrix, u_m_coupler,dt*fctp->theta,fctp->iteration_matrix,dt,sourceN);  
    /* z_m=b-z_m */  
    Paso_Update(n,-1.,z_m,1.,b);  
   
    Paso_Coupler_finishCollect(RN_m_coupler);  
    Paso_Coupler_finishCollect(RP_m_coupler);  
    /* add corrected fluxes into z_m */  
    Paso_FCTransportProblem_addCorrectedFluxes(z_m,flux_matrix,RN_m_coupler,RP_m_coupler);  
    return NO_ERROR;  
 }  
Legend:
Removed from v.1811  
changed lines
  Added in v.1812
  ViewVC Help
Powered by ViewVC 1.1.26