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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 8154 byte(s)
Assorted spelling/comment fixes in paso.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
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
17 /* Paso: Transport solver
18 *
19 * solves Mu_t=Ku+q
20 *
21 * using an operator splitting approach (K=L+D where L is row sum zero and D is a diagonal matrix):
22 *
23 * - Mu_t=Du+q u(0)=u (reactive part)
24 * - Mv_t=Lv v(0)=u(dt/2) (transport part uses flux correction (FCT))
25 * - Mu_t=Du+q u(dt/2)=v(dt/2) (reactive part)
26 *
27 * to return u(dt)
28 *
29 */
30 /**************************************************************/
31
32 /* Author: l.gross@uq.edu.au */
33
34 /**************************************************************/
35
36 #include "Transport.h"
37 #include "FCTSolver.h"
38 #include "ReactiveSolver.h"
39 #include "Solver.h"
40 #include "PasoUtil.h"
41
42
43 void Paso_TransportProblem_solve(Paso_TransportProblem* fctp, double* u, double dt, double* u0, double* q, Paso_Options* options) {
44 const double reduction_after_divergence_factor = 0.5;
45 const dim_t num_failures_max=50;
46 /* const double increase_after_convergence_factor = 1.1; has not much of an impact if substepping is being used */
47
48 Paso_Performance pp;
49 Paso_Function *fctfunction=NULL;
50 Paso_ReactiveSolver *rsolver=NULL;
51
52 dim_t i_substeps, n_substeps, num_failures=0;
53 double *u_save=NULL;
54 double dt_max=LARGE_POSITIVE_FLOAT, dt2,t=0, dt3;
55 err_t errorCode=SOLVER_NO_ERROR;
56 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
57 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);
58 options->time_step_backtracking_used=FALSE;
59 options->num_iter=0;
60
61
62 if (dt<=0.) {
63 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");
64 }
65 if (blockSize>1) {
66 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >1 is not supported.");
67 }
68 if (options->verbose) {
69 if (fctp->useBackwardEuler) {
70 printf("Paso_TransportProblem_solve: Backward Euler is used (dt = %e)\n",dt);
71 } else {
72 printf("Paso_TransportProblem_solve: Crank-Nicolson is used (dt = %e).\n",dt);
73 }
74 }
75 if (Esys_noError()) {
76 dt_max=Paso_TransportProblem_getSafeTimeStepSize(fctp);
77 /*
78 * allocate memory
79 *
80 */
81 fctfunction=Paso_FCTSolver_Function_alloc(fctp,options);
82 rsolver=Paso_ReactiveSolver_alloc(fctp);
83 u_save=TMPMEMALLOC(n,double);
84 Esys_checkPtr(u_save);
85 }
86 if (Esys_noError()) {
87 /*
88 * let the show begin!!!!
89 *
90 */
91
92 /* while(i_substeps<n_substeps && Esys_noError()) */
93 if (dt_max < LARGE_POSITIVE_FLOAT) {
94 dt2=MIN(dt_max,dt);
95 } else {
96 dt2=dt;
97 }
98 num_failures=0;
99
100 Paso_Copy(n,u,u0); /* copy initial value */
101 Paso_Copy(n,u_save, u); /* create copy for restart in case of failure */
102
103
104 while( (dt-t)>dt*sqrt(EPSILON) && Esys_noError()) {
105
106 n_substeps=ceil((dt-t)/dt2);
107 dt3=(dt-t)/n_substeps;
108 if (options->verbose) printf("Paso_TransportProblem_solve: number of substeps = %d with dt = %e.\n",n_substeps,dt3);
109 i_substeps=0;
110 /* initialize the iteration matrix */
111 Paso_FCTSolver_Function_initialize(dt3, fctp,options, &pp);
112 Paso_ReactiveSolver_initialize(dt3, fctp, options);
113 /* start iteration */
114 for (i_substeps =0; (i_substeps<n_substeps) && (errorCode==SOLVER_NO_ERROR) && Esys_noError(); i_substeps++) {
115 if (options->verbose) printf("Paso_TransportProblem_solve: substep %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3);
116 /* updates u */
117 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt3/2,q, options, &pp); /* Mu_t=Du+q u(0)=u */
118 if (errorCode == SOLVER_NO_ERROR) {
119 errorCode=Paso_FCTSolver_solve(fctfunction, u,dt3,options, &pp); /* Mv_t=Lv v(0)=u(dt/2) */
120 }
121 if (errorCode == SOLVER_NO_ERROR) {
122 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt3/2,q, options, &pp); /* Mu_t=Du+q u(dt/2)=v(dt/2) */
123 }
124 if (errorCode == SOLVER_NO_ERROR) {
125 t+=dt3;
126 Paso_Copy(n,u_save, u); /* copy u to u_save */
127 }
128 }
129 /*
130 * error handling:
131 */
132 if (errorCode == SOLVER_NO_ERROR) {
133 num_failures=0;
134 i_substeps++;
135 /*
136 if (fctp->dt_max < LARGE_POSITIVE_FLOAT) {
137 fctp->dt_max=MIN3(dt_max, dt3*increase_after_convergence_factor ,fctp->dt_failed/increase_after_convergence_factor);
138 } else {
139 fctp->dt_max=MIN(dt3*increase_after_convergence_factor ,fctp->dt_failed/increase_after_convergence_factor);
140 }
141 */
142 } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) {
143 /* if num_failures_max failures in a row: give up */
144 fctp->dt_failed=MIN(fctp->dt_failed, dt3);
145 if (num_failures >= num_failures_max) {
146 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: No convergence after time step reductions.");
147 } else {
148 options->time_step_backtracking_used=TRUE;
149 if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n");
150 dt2=dt3*reduction_after_divergence_factor;
151 num_failures++;
152 Paso_Copy(n,u, u_save); /* reset initial value */
153 errorCode = SOLVER_NO_ERROR;
154 }
155
156 } else if (errorCode == SOLVER_INPUT_ERROR) {
157 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver.");
158 } else if (errorCode == SOLVER_MEMORY_ERROR) {
159 Esys_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allocation failed.");
160 } else if (errorCode == SOLVER_BREAKDOWN) {
161 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down.");
162 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
163 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm.");
164 } else {
165 Esys_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error.");
166 }
167 }
168 }
169 /*
170 * clean-up:
171 *
172 */
173 Paso_FCTSolver_Function_free(fctfunction);
174 Paso_ReactiveSolver_free(rsolver);
175 TMPMEMFREE(u_save);
176
177 }
178
179
180 double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp)
181 {
182 double dt_max, dt1=LARGE_POSITIVE_FLOAT, dt2=LARGE_POSITIVE_FLOAT;
183 if ( ! fctp->valid_matrices) {
184 fctp->dt_failed=LARGE_POSITIVE_FLOAT;
185 /* set row-sum of mass_matrix */
186 Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);
187 /* split off row-sum from transport_matrix */
188 Paso_SystemMatrix_makeZeroRowSums(fctp->transport_matrix,fctp->reactive_matrix);
189 /* get a copy of the main diagonal of the mass matrix */
190 Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix);
191
192 if (Esys_noError()) {
193 dt1=Paso_ReactiveSolver_getSafeTimeStepSize(fctp);
194 if (dt1<LARGE_POSITIVE_FLOAT) dt1*=2;
195 }
196 if (Esys_noError()) dt2=Paso_FCTSolver_getSafeTimeStepSize(fctp);
197 printf("Paso_TransportProblem_getSafeTimeStepSize: dt_max from reactive part = %e\n",dt1);
198 printf("Paso_TransportProblem_getSafeTimeStepSize: dt_max from transport part = %e\n",dt2);
199 dt_max=MIN(dt1,dt2);
200
201 if ( (dt_max <= 0.) && Esys_noError() ) {
202 Esys_setError(TYPE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");
203 }
204 fctp->dt_max=dt_max;
205 fctp->valid_matrices=Esys_noError();
206 }
207 return fctp->dt_max;
208 }

  ViewVC Help
Powered by ViewVC 1.1.26