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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1375 - (show annotations)
Wed Jan 9 00:15:05 2008 UTC (11 years, 6 months ago) by gross
File MIME type: text/plain
File size: 6892 byte(s)
bug in interpolation at reduced face elements fixed.
1 /* $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 #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) {
88
89 index_t i, j;
90 int n_substeps,n;
91 double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t;
92 dim_t n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);
93 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 /* add the advective part + artificial diffusion to the diffusive part = transport-matrix */
101 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<n_rows;++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<n_rows;++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 critical
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 }
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<n_rows;++i) {
174 fctp->transport_matrix_diagonal[i]=
175 fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
176 }
177 } else {
178 /* n_substeps=1; */
179 /* 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 for (i=0;i<n_rows;++i) {
187 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
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 t+=dt2;
198 n++;
199 }
200 }
201 if (Paso_noError()) {
202 #pragma omp parallel for schedule(static) private(i)
203 for (i=0;i<n_rows;++i) {
204 u[i]=fctp->u[i];
205 }
206 }
207 }

  ViewVC Help
Powered by ViewVC 1.1.26