/[escript]/trunk/paso/src/Transport_solve.c
ViewVC logotype

Diff of /trunk/paso/src/Transport_solve.c

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

revision 3013 by gross, Thu Apr 22 05:59:31 2010 UTC revision 3014 by gross, Wed Apr 28 04:05:21 2010 UTC
# Line 41  Line 41 
41    
42    
43  void Paso_TransportProblem_solve(Paso_TransportProblem* fctp, double* u, double dt, double* u0, double* q, Paso_Options* options) {  void Paso_TransportProblem_solve(Paso_TransportProblem* fctp, double* u, double dt, double* u0, double* q, Paso_Options* options) {
44       const double reduction_after_divergence_factor = 0.5;
45       const dim_t num_failures_max=50;
46       const double increase_after_convergence_factor = 1.1; /* has not much of an impact if substepping is beeing used */
47      
48     Paso_Performance pp;       Paso_Performance pp;  
49     Paso_Function *fctfunction=NULL;     Paso_Function *fctfunction=NULL;
50     Paso_ReactiveSolver *rsolver=NULL;     Paso_ReactiveSolver *rsolver=NULL;
51     const dim_t FAILURES_MAX=100;  
52     dim_t i_substeps, Failed=0;     dim_t i_substeps, n_substeps, num_failures=0;
53     double *u_save=NULL;     double *u_save=NULL;
54     double dt_max=LARGE_POSITIVE_FLOAT, dt2,t=0;     double dt_max=LARGE_POSITIVE_FLOAT, dt2,t=0, dt3;
55     err_t errorCode=SOLVER_NO_ERROR;     err_t errorCode=SOLVER_NO_ERROR;
56     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
57     const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);     const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);
# Line 59  void Paso_TransportProblem_solve(Paso_Tr Line 63  void Paso_TransportProblem_solve(Paso_Tr
63     if (blockSize>1) {     if (blockSize>1) {
64         Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >0 is not supported.");         Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >0 is not supported.");
65      }      }
66       if (options->verbose) {
67          if (fctp->useBackwardEuler) {
68             printf("Paso_TransportProblem_solve: Backward Euler is used (dt = %e)\n",dt);
69          } else {
70             printf("Paso_TransportProblem_solve: Crank-Nicolson is used (dt = %e).\n",dt);
71          }
72       }
73     if (Paso_noError()) {     if (Paso_noError()) {
74          dt_max=Paso_TransportProblem_getSafeTimeStepSize(fctp);          dt_max=Paso_TransportProblem_getSafeTimeStepSize(fctp);
75          /*          /*
# Line 75  void Paso_TransportProblem_solve(Paso_Tr Line 86  void Paso_TransportProblem_solve(Paso_Tr
86          * let the show begin!!!!          * let the show begin!!!!
87          *          *
88          */          */
         i_substeps=0;  
89            
90          /* while(i_substeps<n_substeps && Paso_noError()) */          /* while(i_substeps<n_substeps && Paso_noError()) */
91          if (dt_max < LARGE_POSITIVE_FLOAT) {          if (dt_max < LARGE_POSITIVE_FLOAT) {
92              dt2=MIN(dt_max,dt);              dt2=MIN(dt_max,dt);
93          } else {          } else {
94               dt2=dt;              dt2=dt;
95          }          }
96      Failed=0;      num_failures=0;
97            
98      Paso_Copy(n,u,u0);      Paso_Copy(n,u,u0); /* copy initial value */
99      Paso_Copy(n,u_save, u);      Paso_Copy(n,u_save, u); /* create copy for restart in case of failure */
100      while( (t<dt*(1.-sqrt(EPSILON))) && Paso_noError()) {      
101        
102          if (options->verbose) printf("substep step %d at t=%e (step size= %e)\n",i_substeps+1,t+dt2,dt2);      while( (dt-t)>dt*sqrt(EPSILON) && Paso_noError()) {
103            
104          /* updates u */          n_substeps=ceil((dt-t)/dt2);
105          errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt2/2,q, options, &pp); /* Mu_t=Du+q u(0)=u */          dt3=(dt-t)/n_substeps;
106          if (errorCode == NO_ERROR) {            if (options->verbose) printf("Paso_TransportProblem_solve: number of substeps = %d with dt = %e.\n",n_substeps,dt3);
107              errorCode=Paso_FCTSolver_solve(fctfunction, u,dt2,options, &pp); /* Mv_t=Lv   v(0)=u(dt/2) */              i_substeps=0;
108          }          /* initialize the iteration matrix */
109          if (errorCode == NO_ERROR) {                Paso_FCTSolver_Function_initialize(dt3, fctp,options, &pp);
110                  errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt2/2,q, options, &pp); /*  Mu_t=Du+q u(dt/2)=v(dt/2) */              Paso_ReactiveSolver_initialize(dt3, fctp, options);
111          }          /* start iteration */
112            for (i_substeps =0; (i_substeps<n_substeps) && (errorCode==SOLVER_NO_ERROR) && Paso_noError(); i_substeps++) {
113                if (options->verbose) printf("Paso_TransportProblem_solve: substep step %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3);
114                /* updates u */
115                errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt3/2,q, options, &pp); /* Mu_t=Du+q u(0)=u */
116                if (errorCode == SOLVER_NO_ERROR) {
117                     errorCode=Paso_FCTSolver_solve(fctfunction, u,dt3,options, &pp); /* Mv_t=Lv   v(0)=u(dt/2) */
118                }
119                if (errorCode == SOLVER_NO_ERROR) {  
120                       errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt3/2,q, options, &pp); /*  Mu_t=Du+q u(dt/2)=v(dt/2) */
121                }
122            if (errorCode == SOLVER_NO_ERROR) {
123                        t+=dt3;
124                        Paso_Copy(n,u_save, u); /* copy u to u_save */
125            }
126                }
127              /*              /*
128           * error handeling:           * error handeling:
129           */           */
130              if (errorCode == SOLVER_NO_ERROR) {               if (errorCode == SOLVER_NO_ERROR) {
131                      Failed=0;                      num_failures=0;
132                      i_substeps++;                      i_substeps++;
                     t+=dt2;  
133                      if (fctp->dt_max < LARGE_POSITIVE_FLOAT) {                      if (fctp->dt_max < LARGE_POSITIVE_FLOAT) {
134                         dt2=MIN3(dt_max,dt2*1.1,dt-t);                         fctp->dt_max=MIN3(dt_max, dt3*increase_after_convergence_factor ,fctp->dt_failed/increase_after_convergence_factor);
135                      } else {                      } else {
136                         dt2=MIN(dt2*1.1,dt-t);                         fctp->dt_max=MIN(dt3*increase_after_convergence_factor ,fctp->dt_failed/increase_after_convergence_factor);
137                      }                      }
138              Paso_Copy(n,u_save, u);           } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) {
139              } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) {                      /* if num_failures_max failures in a row: give up */
140                      /* if FAILURES_MAX failures in a row: give up */              fctp->dt_failed=MIN(fctp->dt_failed, dt3);
141                      if (Failed > FAILURES_MAX) {                      if (num_failures >= num_failures_max) {
142                         Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: no convergence after time step reduction.");                         Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: No convergence after time step reductions.");
143                      } else {                      } else {
144                         if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n");                         if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n");
145                         dt2*=0.5;                         dt2=dt3*reduction_after_divergence_factor;
146                         Failed++;                         num_failures++;
147                 Paso_Copy(n,u, u_save);                 Paso_Copy(n,u, u_save); /* reset initial value */
148                   errorCode = SOLVER_NO_ERROR;
149                      }                      }
150    
151              } else if (errorCode == SOLVER_INPUT_ERROR) {              } else if (errorCode == SOLVER_INPUT_ERROR) {
152              Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver.");              Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver.");
153          } else if (errorCode == SOLVER_MEMORY_ERROR) {          } else if (errorCode == SOLVER_MEMORY_ERROR) {
# Line 144  void Paso_TransportProblem_solve(Paso_Tr Line 170  void Paso_TransportProblem_solve(Paso_Tr
170      TMPMEMFREE(u_save);      TMPMEMFREE(u_save);
171        
172  }  }
173    
174    
175  double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp)  double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp)
176  {  {
177     double dt_max, dt1, dt2;     double dt_max, dt1, dt2;
178     if ( ! fctp->valid_matrices) {     if ( ! fctp->valid_matrices) {
179         fctp->dt_failed=LARGE_POSITIVE_FLOAT;
180       /* set row-sum of mass_matrix */       /* set row-sum of mass_matrix */
181       Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);       Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);
182       /* split off row-sum from transport_matrix */       /* split off row-sum from transport_matrix */

Legend:
Removed from v.3013  
changed lines
  Added in v.3014

  ViewVC Help
Powered by ViewVC 1.1.26