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

Annotation of /trunk/paso/src/SolverFCT_solve.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 6 months ago) by phornby
Original Path: temp_trunk_copy/paso/src/SolverFCT_solve.c
File MIME type: text/plain
File size: 6892 byte(s)
Make a temp copy of the trunk before checking in the windows changes


1 gross 1363 /* $Id:$ */
2    
3     /*******************************************************
4     *
5     * Copyright 2007 by University of Queensland
6     *
7     * http://esscc.uq.edu.au
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14     /**************************************************************/
15    
16     /* Paso: Flux correction transport solver
17     *
18     * solves Mu_t=Du+Ku+q
19     *
20     * where is D is diffusive (not checked)
21     * - D is symmetric
22     * - row sums are equal to zero.
23     * and K is the advective part.
24     *
25     * u(0) >= 0
26     *
27     * intially fctp->transport_matrix defines the diffusive part
28     * but the matrix is updated by the adevctive part + artificial diffusion
29     *
30     */
31     /**************************************************************/
32    
33     /* Author: l.gross@uq.edu.au */
34    
35     /**************************************************************/
36    
37     #include "Paso.h"
38     #include "Solver.h"
39     #include "SolverFCT.h"
40     #include "escript/blocktimer.h"
41 gross 1370 #ifdef _OPENMP
42     #include <omp.h>
43     #endif
44     #ifdef PASO_MPI
45     #include <mpi.h>
46     #endif
47 gross 1363
48 gross 1370 /***********************************************************************************
49 gross 1363
50 gross 1370 solve (until convergence)
51 gross 1363
52 gross 1370 (*) [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 gross 1364 void Paso_SolverFCT_solve(Paso_FCTransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options) {
88 gross 1363
89 gross 1375 index_t i, j;
90 gross 1370 int n_substeps,n;
91     double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t;
92 gross 1372 dim_t n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);
93 gross 1363 if (dt<=0.) {
94     Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
95     }
96     if (! fctp->valid_matrices) {
97    
98     /* extract the row sum of the advective part */
99     Paso_SystemMatrix_rowSum(fctp->flux_matrix,fctp->row_sum_flux_matrix);
100 gross 1370 /* add the advective part + artificial diffusion to the diffusive part = transport-matrix */
101 gross 1363 Paso_FCTransportProblem_addAdvectivePart(fctp,1.);
102 gross 1370 /* create a copy of the main diagonal entires of the transport-matrix */
103     #pragma omp parallel for schedule(static) private(i)
104 gross 1372 for (i=0;i<n_rows;++i) {
105 gross 1370 fctp->transport_matrix_diagonal[i]=
106     fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
107     }
108 gross 1363
109 gross 1370 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 gross 1372 for (i=0;i<n_rows;++i) {
123 gross 1370 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 gross 1372 #pragma omp critical
130 gross 1370 {
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     }
158     fctp->valid_matrices=Paso_noError();
159 gross 1363 }
160    
161 gross 1370 /* 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 gross 1372 for (i=0;i<n_rows;++i) {
174 gross 1370 fctp->transport_matrix_diagonal[i]=
175     fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
176     }
177     } else {
178 gross 1371 /* n_substeps=1; */
179 gross 1370 /* u= u_last + M^-1*dt*F(u_last) */
180     t=dt2;
181     n=0;
182     while(n<n_substeps) {
183     printf("substep step %d at t=%e\n",n+1,t);
184     Paso_FCTransportProblem_setFlux(fctp,fctp->u,u); /* u stores F(u_last)*/
185     #pragma omp parallel for schedule(static) private(i)
186 gross 1372 for (i=0;i<n_rows;++i) {
187 gross 1370 rtmp=fctp->u[i];
188     fctp->u[i]=fctp->u[i]+dt2*u[i]/fctp->lumped_mass_matrix[i];
189     printf("%d : %e %e %e\n",i,u[i],fctp->u[i],rtmp);
190     }
191 gross 1375
192     for (i=0;i<21;++i) {
193     for (j=0;j<21;++j) printf("%d->%e ",i*21+j,fctp->u[i*21+j]);
194     printf("\n");
195     }
196    
197 gross 1370 t+=dt2;
198     n++;
199     }
200     }
201     if (Paso_noError()) {
202     #pragma omp parallel for schedule(static) private(i)
203 gross 1372 for (i=0;i<n_rows;++i) {
204 gross 1370 u[i]=fctp->u[i];
205     }
206     }
207     }

  ViewVC Help
Powered by ViewVC 1.1.26