/[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 1371 - (show annotations)
Thu Jan 3 06:11:21 2008 UTC (11 years, 7 months ago) by gross
File MIME type: text/plain
File size: 6941 byte(s)
some bugs in the updwing scheme but there are still problem with the cross wind direction.
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;
90 int n_substeps,n;
91 double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t;
92 if (dt<=0.) {
93 Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive.");
94 }
95 if (! fctp->valid_matrices) {
96
97 /* extract the row sum of the advective part */
98 Paso_SystemMatrix_rowSum(fctp->flux_matrix,fctp->row_sum_flux_matrix);
99 /* add the advective part + artificial diffusion to the diffusive part = transport-matrix */
100 Paso_FCTransportProblem_addAdvectivePart(fctp,1.);
101 /* create a copy of the main diagonal entires of the transport-matrix */
102 #pragma omp parallel for schedule(static) private(i)
103 for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {
104 fctp->transport_matrix_diagonal[i]=
105 fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
106 }
107
108 Paso_SystemMatrix_saveMM(fctp->flux_matrix,"flux.mm");
109 Paso_SystemMatrix_saveMM(fctp->transport_matrix,"trans.mm");
110 /*=================================================================== *
111 *
112 * calculate time step size:
113 */
114 dt2=fctp->dt_max;
115 if (fctp->theta < 1.) {
116 dt2=LARGE_POSITIVE_FLOAT;
117 #pragma omp parallel private(dt2_loc)
118 {
119 dt2_loc=LARGE_POSITIVE_FLOAT;
120 #pragma omp for schedule(static) private(i,rtmp,rtmp2)
121 for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {
122 rtmp=fctp->transport_matrix_diagonal[i];
123 rtmp2=fctp->lumped_mass_matrix[i];
124 if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) {
125 dt2_loc=MIN(dt2_loc,-rtmp2/rtmp);
126 }
127 }
128 #pragma omp critcal
129 {
130 dt2=MIN(dt2,dt2_loc);
131 }
132 }
133 #ifdef PASO_MPI
134 dt2_loc = dt2;
135 MPI_Allreduce(&dt2_loc, &dt2, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
136 #endif
137 dt2*=1./(1.-fctp->theta);
138 if (fctp->dt_max>0.) dt2=MIN(dt2,fctp->dt_max);
139 }
140 if (dt2 > 0.) {
141 fctp->dt=dt2;
142 } else {
143 fctp->dt=fctp->dt_max;
144 }
145 if (fctp->dt < 0.) {
146 Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt_max must be positive.");
147 } else {
148 printf("minimum time step size is %e (theta = %e).\n",fctp->dt,fctp->theta);
149 }
150 /* ===========================================================
151 *
152 *
153 */
154 if (Paso_noError()) {
155
156 }
157 fctp->valid_matrices=Paso_noError();
158 }
159
160 /* b_last=M*u_last + (1-theta) * F(u_last) */
161
162 /* decide on substepping */
163 n_substeps=ceil(dt/fctp->dt);
164 dt2=dt/n_substeps;
165 printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta,fctp->dt);
166
167 /* */
168 Paso_SystemMatrix_allocBuffer(fctp->transport_matrix);
169
170 if (fctp->theta>0) {
171 #pragma omp parallel for schedule(static) private(i)
172 for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) {
173 fctp->transport_matrix_diagonal[i]=
174 fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]];
175 }
176 } else {
177 /* n_substeps=1; */
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 }

  ViewVC Help
Powered by ViewVC 1.1.26