/[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 1369 by gross, Mon Dec 17 07:22:45 2007 UTC revision 1370 by gross, Wed Jan 2 09:21:43 2008 UTC
# Line 38  Line 38 
38  #include "Solver.h"  #include "Solver.h"
39  #include "SolverFCT.h"  #include "SolverFCT.h"
40  #include "escript/blocktimer.h"  #include "escript/blocktimer.h"
41    #ifdef _OPENMP
42    #include <omp.h>
43    #endif
44    #ifdef PASO_MPI
45    #include <mpi.h>
46    #endif
47    
48    /***********************************************************************************
49    
50            solve (until convergence)
51    
52            (*)  [M-theta*dt*L] du = M*u_last + (1-theta)*dt*F(u_last) - M*u + theta*dt F(u)
53                 u <- u+du
54    
55            with F(u) =  L u + f_a(u)  (f_a anti-diffusion)
56            and L = D + K + D_K  stored in transport_matrix
57                D = Diffusive part (on input stored in transport_matrix)
58                K = flux_matrix
59                D_K = artificial diffusion introduced by K
60                f_a = anti-diffusion introduced by K and u
61                u_last=u of last time step
62    
63            For the case theta==0 (explicit case) no iteration is required. One sets
64    
65                  M*u= M*u_last + dt*F(u_last)
66    
67            with F(u) =  L u + f_a(u).
68    
69    
70            For the case theta>0 we set
71    
72                A=L-M/(theta*dt) (stored into transport_matrix)
73                F(u)=L u + f_a(u) = M/(theta*dt) u + A u + f_a(u)
74                b(u)= - M/(theta*dt) * u + F(u) =  A u + f_a(u)
75                b_last=M/(theta*dt)*u_last + (1-theta)*dt * F(u_last)
76                      =M/(theta*dt)*u_last + (1-theta)/theta * [ M/(theta*dt) u_last + A u_last + f_a(u_last) ]
77                      =(1/theta**2*dt)*M u_last + (1-theta)/theta*b(u_last)
78            so (*) takes the form
79    
80            b_last=(1/theta**2*dt)*M u_last + (1-theta)/theta*b(u_last)
81            while (|du| > tol * |u|) :
82                 A*du=b_last + b(u)
83                 u <- u-du      
84    
85  /***********************************************************************************/  */
   
86    
87  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) {
88    
89       index_t i;
90       int n_substeps,n;
91       double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t;
92       printf("FFFF %e",fctp->dt_max);
93     if (dt<=0.) {     if (dt<=0.) {
94         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");         Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
95     }     }
# Line 51  void Paso_SolverFCT_solve(Paso_FCTranspo Line 97  void Paso_SolverFCT_solve(Paso_FCTranspo
97    
98          /* extract the row sum of the advective part */          /* extract the row sum of the advective part */
99          Paso_SystemMatrix_rowSum(fctp->flux_matrix,fctp->row_sum_flux_matrix);          Paso_SystemMatrix_rowSum(fctp->flux_matrix,fctp->row_sum_flux_matrix);
100            /* add the advective part + artificial diffusion to the diffusive part = transport-matrix */
         /* add the advective part + artificial diffusion to the diffusive part */  
101          Paso_FCTransportProblem_addAdvectivePart(fctp,1.);          Paso_FCTransportProblem_addAdvectivePart(fctp,1.);
102            /* create a copy of the main diagonal entires of the transport-matrix */
103            #pragma omp parallel for schedule(static) private(i)
104            for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {
105                   fctp->transport_matrix_diagonal[i]=
106                            fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
107            }
108    
109    Paso_SystemMatrix_saveMM(fctp->flux_matrix,"flux.mm");
110    Paso_SystemMatrix_saveMM(fctp->transport_matrix,"trans.mm");
111            /*=================================================================== *
112             *
113             *  calculate time step size:                                          
114             */
115            dt2=fctp->dt_max;
116            if (fctp->theta < 1.) {
117                dt2=LARGE_POSITIVE_FLOAT;
118                #pragma omp parallel private(dt2_loc)
119                {
120                   dt2_loc=LARGE_POSITIVE_FLOAT;
121                   #pragma omp for schedule(static) private(i,rtmp,rtmp2)
122                   for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {
123                        rtmp=fctp->transport_matrix_diagonal[i];
124                        rtmp2=fctp->lumped_mass_matrix[i];
125                        if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) {
126                            dt2_loc=MIN(dt2_loc,-rtmp2/rtmp);
127                        }
128                    }
129                    #pragma omp critcal
130                    {
131                        dt2=MIN(dt2,dt2_loc);
132                    }
133                }
134                #ifdef PASO_MPI
135                   dt2_loc = dt2;
136               MPI_Allreduce(&dt2_loc, &dt2, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
137                #endif
138                dt2*=1./(1.-fctp->theta);
139                if (fctp->dt_max>0.) dt2=MIN(dt2,fctp->dt_max);
140            }
141            if (dt2 > 0.) {
142             fctp->dt=dt2;
143            } else {
144             fctp->dt=fctp->dt_max;
145            }
146            if (fctp->dt < 0.) {
147                Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt_max must be positive.");
148            } else {
149               printf("minimum time step size is %e (theta = %e).\n",fctp->dt,fctp->theta);  
150            }
151            /* ===========================================================
152             *
153             *    
154            */
155            if (Paso_noError()) {
156    
157          if (Paso_noError()) fctp->valid_matrices=TRUE;          }
158            fctp->valid_matrices=Paso_noError();
159     }     }
160    
161  }     /* b_last=M*u_last + (1-theta) * F(u_last) */
162    
163       /* decide on substepping */
164       n_substeps=ceil(dt/fctp->dt);
165       dt2=dt/n_substeps;
166       printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta,fctp->dt);
167    
168       /* */
169       Paso_SystemMatrix_allocBuffer(fctp->transport_matrix);
170    
171       if (fctp->theta>0) {
172            #pragma omp parallel for schedule(static) private(i)
173            for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {
174                   fctp->transport_matrix_diagonal[i]=
175                            fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
176            }
177       } else {
178          /* u= u_last + M^-1*dt*F(u_last) */
179          t=dt2;
180          n=0;
181          while(n<n_substeps) {
182            printf("substep step %d at t=%e\n",n+1,t);
183            Paso_FCTransportProblem_setFlux(fctp,fctp->u,u); /* u stores F(u_last)*/
184            #pragma omp parallel for schedule(static) private(i)
185            for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {
186                rtmp=fctp->u[i];
187                fctp->u[i]=fctp->u[i]+dt2*u[i]/fctp->lumped_mass_matrix[i];
188                printf("%d : %e %e %e\n",i,u[i],fctp->u[i],rtmp);
189            }
190            t+=dt2;
191            n++;
192          }
193      }
194      if (Paso_noError()) {
195         #pragma omp parallel for schedule(static) private(i)
196         for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {
197           u[i]=fctp->u[i];
198         }
199      }
200    }      
Legend:
Removed from v.1369  
changed lines
  Added in v.1370
  ViewVC Help
Powered by ViewVC 1.1.26