/[escript]/branches/stage3.1/paso/src/SolverFCT_solve.c
ViewVC logotype

Diff of /branches/stage3.1/paso/src/SolverFCT_solve.c

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

revision 2944 by jfenwick, Thu Feb 4 01:42:47 2010 UTC revision 2945 by jfenwick, Wed Feb 24 00:17:46 2010 UTC
# Line 98  err_t Paso_FCT_setUpRightHandSide(Paso_F Line 98  err_t Paso_FCT_setUpRightHandSide(Paso_F
98       */       */
99     Paso_SolverFCT_setMuPaLuPbQ(z_m,fctp->lumped_mass_matrix, u_m_coupler,dt*fctp->theta,fctp->iteration_matrix,dt,sourceN);     Paso_SolverFCT_setMuPaLuPbQ(z_m,fctp->lumped_mass_matrix, u_m_coupler,dt*fctp->theta,fctp->iteration_matrix,dt,sourceN);
100     /* z_m=b-z_m */     /* z_m=b-z_m */
101     /*{
102     int kk;
103     for (kk=0;kk<n;kk++) printf("z_m %d : %e %e -> %e\n",kk,z_m[kk], b[kk],z_m[kk]-b[kk]);
104     }
105    */
106     Paso_Update(n,-1.,z_m,1.,b);     Paso_Update(n,-1.,z_m,1.,b);
107    
108     Paso_Coupler_finishCollect(RN_m_coupler);     Paso_Coupler_finishCollect(RN_m_coupler);
109     Paso_Coupler_finishCollect(RP_m_coupler);     Paso_Coupler_finishCollect(RP_m_coupler);
110     /* add corrected fluxes into z_m */     /* add corrected fluxes into z_m */
111     Paso_FCTransportProblem_addCorrectedFluxes(z_m,flux_matrix,RN_m_coupler,RP_m_coupler);      Paso_FCTransportProblem_addCorrectedFluxes(z_m,flux_matrix,RN_m_coupler,RP_m_coupler);
112     return NO_ERROR;     return NO_ERROR;
113  }  }
114    
# Line 198  void Paso_SolverFCT_solve(Paso_FCTranspo Line 203  void Paso_SolverFCT_solve(Paso_FCTranspo
203                       * now we solve the linear system to get the correction dt:                       * now we solve the linear system to get the correction dt:
204                       *                       *
205                       */                       */
 /* for (i=0;i<n;++i) printf("%d %f \n", i,z_m[i]); */  
206                       if (fctp->theta > 0) {                       if (fctp->theta > 0) {
207                            omega=1./(dt2*fctp->theta);                            omega=1./(dt2*fctp->theta);
208                            Paso_Solver_solvePreconditioner(fctp->iteration_matrix,du_m,z_m);                              Paso_Solver_solvePreconditioner(fctp->iteration_matrix,du_m,z_m);  
209                            Paso_Update(n,1.,u,omega,du_m);                            Paso_Update(n,1.,u,omega,du_m);
210                       } else {  /* for (i = 0; i < n; ++i) printf("u %d = %e %e <= %e \n",i,u[i], du_m[i], z_m[i]); */
211                 } else {
212                            omega=1;                            omega=1;
213                            #pragma omp parallel for private(i,mass,rtmp)                            #pragma omp parallel for private(i,mass,rtmp)
214                            for (i = 0; i < n; ++i) {                            for (i = 0; i < n; ++i) {
# Line 216  void Paso_SolverFCT_solve(Paso_FCTranspo Line 221  void Paso_SolverFCT_solve(Paso_FCTranspo
221                                du_m[i]=rtmp;                                du_m[i]=rtmp;
222                                u[i]+=rtmp;                                u[i]+=rtmp;
223                            }                            }
224                     }                       }
225                     norm_u_m=Paso_lsup(n,u, fctp->mpi_info);                       norm_u_m=Paso_lsup(n,u, fctp->mpi_info);
226                     norm_du_m=Paso_lsup(n,du_m, fctp->mpi_info)*omega;                       norm_du_m=Paso_lsup(n,du_m, fctp->mpi_info)*omega;
227                     if (options->verbose) printf("iteration step %d completed: norm increment= %e (tolerance = %e)\n",m+1, norm_du_m, rtol * norm_u_m + atol);                       if (options->verbose) printf("iteration step %d completed: norm increment= %e (tolerance = %e)\n",m+1, norm_du_m, rtol * norm_u_m + atol);
228    
229                     max_m_reached=(m>max_m);                       max_m_reached=(m>max_m);
230                     converged=(norm_du_m <= rtol * norm_u_m + atol);                       converged=(norm_du_m <= rtol * norm_u_m + atol);
231                     m++;                       m++;
232              }              }
233              if (converged) {              if (converged) {
234                      Failed=0;                      Failed=0;
# Line 362  void Paso_FCT_setUp(Paso_FCTransportProb Line 367  void Paso_FCT_setUp(Paso_FCTransportProb
367                if (ABS(u_tilde_i)>0) rtmp4+=sourceN[i]*factor/u_tilde_i;                if (ABS(u_tilde_i)>0) rtmp4+=sourceN[i]*factor/u_tilde_i;
368                fctp->iteration_matrix->mainBlock->val[fctp->main_iptr[i]]=rtmp4;                fctp->iteration_matrix->mainBlock->val[fctp->main_iptr[i]]=rtmp4;
369                uTilde[i]=u_tilde_i;                uTilde[i]=u_tilde_i;
370    /* printf("uTilde %d %e and %e : %e : %e %e\n",i,u_tilde_i, b[i], rtmp4, m, fctp->main_diagonal_low_order_transport_matrix[i]); */
371           }           }
372           Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);           Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
373           Paso_Solver_setPreconditioner(fctp->iteration_matrix,options);           Paso_Solver_setPreconditioner(fctp->iteration_matrix,options);

Legend:
Removed from v.2944  
changed lines
  Added in v.2945

  ViewVC Help
Powered by ViewVC 1.1.26