/[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 3792 by caltinay, Thu Oct 27 03:41:51 2011 UTC revision 3793 by gross, Wed Feb 1 07:39:43 2012 UTC
# Line 33  Line 33 
33    
34  /**************************************************************/  /**************************************************************/
35    
36    
37    
38  #include "Transport.h"  #include "Transport.h"
39  #include "FCTSolver.h"  #include "FCT_Solver.h"
40  #include "ReactiveSolver.h"  #include "ReactiveSolver.h"
41  #include "Solver.h"  #include "Solver.h"
42  #include "PasoUtil.h"  #include "PasoUtil.h"
# Line 43  Line 45 
45  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) {
46     const double reduction_after_divergence_factor = 0.5;     const double reduction_after_divergence_factor = 0.5;
47     const dim_t num_failures_max=50;     const dim_t num_failures_max=50;
    /* const double increase_after_convergence_factor = 1.1;  has not much of an impact if substepping is being used */  
48        
49     Paso_Performance pp;       Paso_Performance pp;
    Paso_Function *fctfunction=NULL;  
50     Paso_ReactiveSolver *rsolver=NULL;     Paso_ReactiveSolver *rsolver=NULL;
51       Paso_FCT_Solver *fctsolver=NULL;
52          
53    
54     dim_t i_substeps, n_substeps, num_failures=0;     dim_t i_substeps, n_substeps, num_failures=0;
55     double *u_save=NULL;     double *u_save=NULL, *u2=NULL;
56     double dt_max=LARGE_POSITIVE_FLOAT, dt2,t=0, dt3;     double  dt2,t=0, dt3;
57     err_t errorCode=SOLVER_NO_ERROR;     err_t errorCode=SOLVER_NO_ERROR;
58     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
59     const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);     const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);
60     options->time_step_backtracking_used=FALSE;     options->time_step_backtracking_used=FALSE;
61     options->num_iter=0;     options->num_iter=0;  
   
62    
63     if (dt<=0.) {     if (dt<=0.) {
64         Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");         Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");
# Line 66  void Paso_TransportProblem_solve(Paso_Tr Line 67  void Paso_TransportProblem_solve(Paso_Tr
67         Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >1 is not supported.");         Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >1 is not supported.");
68      }      }
69     if (options->verbose) {     if (options->verbose) {
70        if (fctp->useBackwardEuler) {        if (options->ode_solver == PASO_BACKWARD_EULER) {
71           printf("Paso_TransportProblem_solve: Backward Euler is used (dt = %e)\n",dt);           printf("Paso_TransportProblem_solve: Backward Euler is used (dt = %e)\n",dt);
72        } else {        } else  if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) {
73             printf("Paso_TransportProblem_solve: linear Crank-Nicolson is used (dt = %e).\n",dt);
74          } else  if (options->ode_solver == PASO_CRANK_NICOLSON) {
75           printf("Paso_TransportProblem_solve: Crank-Nicolson is used (dt = %e).\n",dt);           printf("Paso_TransportProblem_solve: Crank-Nicolson is used (dt = %e).\n",dt);
76          } else {
77         Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: unknown ODE solver.");
78        }        }
79     }     }
80     if (Esys_noError()) {     if (Esys_noError()) {
81          dt_max=Paso_TransportProblem_getSafeTimeStepSize(fctp);           Paso_TransportProblem_getSafeTimeStepSize(fctp);
82          /*          /*
83           *  allocate memory           *  allocate memory
84           *           *
85           */           */
86            fctfunction=Paso_FCTSolver_Function_alloc(fctp,options);            fctsolver=Paso_FCT_Solver_alloc(fctp,options);
87            rsolver=Paso_ReactiveSolver_alloc(fctp);            rsolver=Paso_ReactiveSolver_alloc(fctp);
88        u_save=TMPMEMALLOC(n,double);        u_save=TMPMEMALLOC(n,double);
89          u2=TMPMEMALLOC(n,double);
90            Esys_checkPtr(u_save);            Esys_checkPtr(u_save);
91          Esys_checkPtr(u2);
92      }      }
93     if (Esys_noError()) {       if (Esys_noError()) {  
94         /*         /*
95          * let the show begin!!!!          * let the show begin!!!!
96          *          *
97          */          */
98        const double dt_R=fctp->dt_max_R;
99            const double dt_T=fctp->dt_max_T;
100        dt2=dt;
101        if (dt_R < LARGE_POSITIVE_FLOAT) dt2 = MIN(dt_R *2, dt); /* as we half the steps size for the RT bit */
102        if (dt_T < LARGE_POSITIVE_FLOAT) {
103              if ( ( options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) || ( options->ode_solver == PASO_CRANK_NICOLSON ) ) {
104                 dt2 = MIN(dt_T, dt);
105              } /* PASO_BACKWARD_EULER does not require a restriction */
106        }
107            
108          /* while(i_substeps<n_substeps && Esys_noError()) */          num_failures=0;
109          if (dt_max < LARGE_POSITIVE_FLOAT) {      Paso_Copy(n,u,u0); /* copy initial value to return */
110              dt2=MIN(dt_max,dt);          
         } else {  
             dt2=dt;  
         }  
     num_failures=0;  
       
     Paso_Copy(n,u,u0); /* copy initial value */  
     Paso_Copy(n,u_save, u); /* create copy for restart in case of failure */  
       
       
111      while( (dt-t)>dt*sqrt(EPSILON) && Esys_noError()) {      while( (dt-t)>dt*sqrt(EPSILON) && Esys_noError()) {
112    
113          n_substeps=ceil((dt-t)/dt2);          n_substeps=ceil((dt-t)/dt2);
# Line 108  void Paso_TransportProblem_solve(Paso_Tr Line 115  void Paso_TransportProblem_solve(Paso_Tr
115          if (options->verbose) printf("Paso_TransportProblem_solve: number of substeps = %d with dt = %e.\n",n_substeps,dt3);          if (options->verbose) printf("Paso_TransportProblem_solve: number of substeps = %d with dt = %e.\n",n_substeps,dt3);
116              i_substeps=0;              i_substeps=0;
117          /* initialize the iteration matrix */          /* initialize the iteration matrix */
118              Paso_FCTSolver_Function_initialize(dt3, fctp,options, &pp);          Paso_FCT_Solver_initialize(dt3, fctsolver, options, &pp);
119              Paso_ReactiveSolver_initialize(dt3, fctp, options);              Paso_ReactiveSolver_initialize(dt3/2, rsolver, options);
120            errorCode = SOLVER_NO_ERROR;
121    
122          /* start iteration */          /* start iteration */
123          for (i_substeps =0; (i_substeps<n_substeps) && (errorCode==SOLVER_NO_ERROR) && Esys_noError(); i_substeps++) {          for (i_substeps =0; (i_substeps<n_substeps) && (errorCode==SOLVER_NO_ERROR) && Esys_noError(); i_substeps++) {
124              if (options->verbose) printf("Paso_TransportProblem_solve: substep %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3);              if (options->verbose) printf("Paso_TransportProblem_solve: substep %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3);
125              /* updates u */              Paso_Copy(n,u_save, u); /* create copy for restart in case of failure */
126              errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt3/2,q, options, &pp); /* Mu_t=Du+q u(0)=u */          /* updates u */
127                errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u2, u ,q, options, &pp); /* Mu_t=Du+q u(0)=u */
128              if (errorCode == SOLVER_NO_ERROR) {              if (errorCode == SOLVER_NO_ERROR) {
129                   errorCode=Paso_FCTSolver_solve(fctfunction, u,dt3,options, &pp); /* Mv_t=Lv   v(0)=u(dt/2) */                    errorCode=Paso_FCT_Solver_update(fctsolver, u, u2,options, &pp); /* Mv_t=Lv   v(0)=u(dt/2) */
130                    
131              }              }
132              if (errorCode == SOLVER_NO_ERROR) {                if (errorCode == SOLVER_NO_ERROR) {  
133                     errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt3/2,q, options, &pp); /*  Mu_t=Du+q u(dt/2)=v(dt/2) */                     errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u2, u, q, options, &pp); /*  Mu_t=Du+q u(dt/2)=v(dt/2) */
134              }              }
135          if (errorCode == SOLVER_NO_ERROR) {                  
136                      t+=dt3;  
137                      Paso_Copy(n,u_save, u); /* copy u to u_save */                  /*
138          }               * error handling:
139              }               */
140              /*                  if (errorCode == SOLVER_NO_ERROR) {
          * error handling:  
          */  
              if (errorCode == SOLVER_NO_ERROR) {  
141                      num_failures=0;                      num_failures=0;
142                      i_substeps++;              t+=dt3;
143  /*              Paso_Copy(n,u, u2);
144              if (fctp->dt_max < LARGE_POSITIVE_FLOAT) {          }
145                         fctp->dt_max=MIN3(dt_max, dt3*increase_after_convergence_factor ,fctp->dt_failed/increase_after_convergence_factor);          }
146                      } else {          if (errorCode == SOLVER_NO_ERROR) {
147                         fctp->dt_max=MIN(dt3*increase_after_convergence_factor ,fctp->dt_failed/increase_after_convergence_factor);          } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) {
                     }  
  */  
          } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) {  
148                      /* if num_failures_max failures in a row: give up */                      /* if num_failures_max failures in a row: give up */
             fctp->dt_failed=MIN(fctp->dt_failed, dt3);  
149                      if (num_failures >= num_failures_max) {                      if (num_failures >= num_failures_max) {
150                         Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: No convergence after time step reductions.");                         Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: No convergence after time step reductions.");
151                      } else {                      } else {
# Line 150  void Paso_TransportProblem_solve(Paso_Tr Line 154  void Paso_TransportProblem_solve(Paso_Tr
154                         dt2=dt3*reduction_after_divergence_factor;                         dt2=dt3*reduction_after_divergence_factor;
155                         num_failures++;                         num_failures++;
156                 Paso_Copy(n,u, u_save); /* reset initial value */                 Paso_Copy(n,u, u_save); /* reset initial value */
                errorCode = SOLVER_NO_ERROR;  
157                      }                      }
158                 } else if (errorCode == SOLVER_INPUT_ERROR) {
             } else if (errorCode == SOLVER_INPUT_ERROR) {  
159              Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver.");              Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver.");
160          } else if (errorCode == SOLVER_MEMORY_ERROR) {           } else if (errorCode == SOLVER_MEMORY_ERROR) {
161              Esys_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allocation failed.");              Esys_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allocation failed.");
162          } else if (errorCode == SOLVER_BREAKDOWN) {           } else if (errorCode == SOLVER_BREAKDOWN) {
163              Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down.");              Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down.");
164          } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {           } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
165              Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm.");              Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm.");
166          } else {           } else {
167              Esys_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error.");              Esys_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error.");
168          }           }
169      }      }
170       }       }
171      /*      /*
172       *  clean-up:       *  clean-up:
173       *       *
174       */       */
175      Paso_FCTSolver_Function_free(fctfunction);  
176        Paso_FCT_Solver_free(fctsolver);
177      Paso_ReactiveSolver_free(rsolver);      Paso_ReactiveSolver_free(rsolver);
178      TMPMEMFREE(u_save);      TMPMEMFREE(u_save);
179          TMPMEMFREE(u2);
180  }  }
181    
182    
183  double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp)  double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp)
184  {  {
185     double dt_max, dt1=LARGE_POSITIVE_FLOAT, dt2=LARGE_POSITIVE_FLOAT;     double dt_max=0., dt_R=LARGE_POSITIVE_FLOAT, dt_T=LARGE_POSITIVE_FLOAT;
186     if ( ! fctp->valid_matrices) {     if ( ! fctp->valid_matrices) {
      fctp->dt_failed=LARGE_POSITIVE_FLOAT;  
187       /* set row-sum of mass_matrix */       /* set row-sum of mass_matrix */
188       Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);       Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);
189       /* split off row-sum from transport_matrix */       /* split off row-sum from transport_matrix */
# Line 190  double Paso_TransportProblem_getSafeTime Line 192  double Paso_TransportProblem_getSafeTime
192       Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix);       Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix);
193    
194       if (Esys_noError()) {       if (Esys_noError()) {
195            dt1=Paso_ReactiveSolver_getSafeTimeStepSize(fctp);            dt_R=Paso_ReactiveSolver_getSafeTimeStepSize(fctp);
196            if (dt1<LARGE_POSITIVE_FLOAT) dt1*=2;        dt_T=Paso_FCT_Solver_getSafeTimeStepSize(fctp);
197       }       }
198       if (Esys_noError()) dt2=Paso_FCTSolver_getSafeTimeStepSize(fctp);       if (Esys_noError()) {
199       printf("Paso_TransportProblem_getSafeTimeStepSize: dt_max from reactive part = %e\n",dt1);            fctp->dt_max_R=dt_R;
200       printf("Paso_TransportProblem_getSafeTimeStepSize: dt_max from transport part = %e\n",dt2);            fctp->dt_max_T=dt_T;
201       dt_max=MIN(dt1,dt2);            fctp->valid_matrices=TRUE;
202          dt_max=MIN(dt_R,dt_T);
203       if ( (dt_max <= 0.) &&  Esys_noError() ) {         }
204              Esys_setError(TYPE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");     } else {
205        }          dt_max=MIN(fctp->dt_max_R,fctp->dt_max_T);  
      fctp->dt_max=dt_max;  
      fctp->valid_matrices=Esys_noError();  
206     }     }
207     return fctp->dt_max;     return dt_max;
208  }  }

Legend:
Removed from v.3792  
changed lines
  Added in v.3793

  ViewVC Help
Powered by ViewVC 1.1.26