/[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 3825 - (show annotations)
Tue Feb 14 06:11:58 2012 UTC (7 years, 2 months ago) by gross
File MIME type: text/plain
File size: 9151 byte(s)
don't call MPI withing parallel regions
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
37
38 #include "Transport.h"
39 #include "FCT_Solver.h"
40 #include "ReactiveSolver.h"
41 #include "Solver.h"
42 #include "PasoUtil.h"
43
44
45 void Paso_TransportProblem_solve(Paso_TransportProblem* fctp, double* u, double dt, double* u0, double* q, Paso_Options* options) {
46 const double reduction_after_divergence_factor = 0.5;
47 const dim_t num_failures_max=50;
48
49 Paso_Performance pp;
50 Paso_ReactiveSolver *rsolver=NULL;
51 Paso_FCT_Solver *fctsolver=NULL;
52
53
54 dim_t i_substeps=0, n_substeps=1, num_failures=0;
55 double *u_save=NULL, *u2=NULL;
56 double dt2,t=0, dt3;
57 err_t errorCode=SOLVER_NO_ERROR;
58 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
59 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);
60 options->time_step_backtracking_used=FALSE;
61 options->num_iter=0;
62
63 if (dt<=0.) {
64 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");
65 }
66 if (blockSize>1) {
67 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >1 is not supported.");
68 }
69 if (options->verbose) {
70 if (options->ode_solver == PASO_BACKWARD_EULER) {
71 printf("Paso_TransportProblem_solve: Backward Euler is used (dt = %e)\n",dt);
72 } else if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) {
73 printf("Paso_TransportProblem_solve: linear Crank-Nicolson is used (dt = %e).\n",dt);
74 } else if (options->ode_solver == PASO_CRANK_NICOLSON) {
75 printf("Paso_TransportProblem_solve: Crank-Nicolson is used (dt = %e).\n",dt);
76 } else {
77 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: unknown ODE solver.");
78 }
79 }
80 if (Esys_noError()) {
81 Paso_TransportProblem_getSafeTimeStepSize(fctp);
82 /*
83 * allocate memory
84 *
85 */
86 fctsolver=Paso_FCT_Solver_alloc(fctp,options);
87 rsolver=Paso_ReactiveSolver_alloc(fctp);
88 u_save=TMPMEMALLOC(n,double);
89 u2=TMPMEMALLOC(n,double);
90 Esys_checkPtr(u_save);
91 Esys_checkPtr(u2);
92 }
93 if (Esys_noError()) {
94 /*
95 * let the show begin!!!!
96 *
97 */
98 const double dt_R=fctp->dt_max_R;
99 const double dt_T=fctp->dt_max_T;
100 dt2=dt;
101 if (dt_R < LARGE_POSITIVE_FLOAT) dt2 = MIN(dt_R *2, dt); /* as we half the steps size for the RT bit */
102 if (dt_T < LARGE_POSITIVE_FLOAT) {
103 if ( ( options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) || ( options->ode_solver == PASO_CRANK_NICOLSON ) ) {
104 dt2 = MIN(dt_T, dt);
105 } /* PASO_BACKWARD_EULER does not require a restriction */
106 }
107
108 num_failures=0;
109 Paso_Copy(n,u,u0); /* copy initial value to return */
110
111 while( (dt-t)>dt*sqrt(EPSILON) && Esys_noError()) {
112
113 n_substeps=ceil((dt-t)/dt2);
114 if ( n_substeps <= 0) {
115 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: time stepping break down.");
116 } else {
117 dt3=(dt-t)/n_substeps;
118 if (options->verbose) printf("Paso_TransportProblem_solve: number of substeps = %d with dt = %e.\n",n_substeps,dt3);
119 /* initialize the iteration matrix */
120 Paso_FCT_Solver_initialize(dt3, fctsolver, options, &pp);
121 Paso_ReactiveSolver_initialize(dt3/2, rsolver, options);
122 errorCode = SOLVER_NO_ERROR;
123
124 /* start iteration */
125 for (i_substeps =0; (i_substeps<n_substeps) && (errorCode==SOLVER_NO_ERROR) && Esys_noError(); i_substeps++) {
126 if (options->verbose) printf("Paso_TransportProblem_solve: substep %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3);
127 Paso_Copy(n,u_save, u); /* create copy for restart in case of failure */
128 /* updates u */
129 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u2, u ,q, options, &pp); /* Mu_t=Du+q u(0)=u */
130 if (errorCode == SOLVER_NO_ERROR) {
131 errorCode=Paso_FCT_Solver_update(fctsolver, u, u2,options, &pp); /* Mv_t=Lv v(0)=u(dt/2) */
132
133 }
134 if (errorCode == SOLVER_NO_ERROR) {
135 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u2, u, q, options, &pp); /* Mu_t=Du+q u(dt/2)=v(dt/2) */
136 }
137
138
139 /*
140 * error handling:
141 */
142 if (errorCode == SOLVER_NO_ERROR) {
143 num_failures=0;
144 t+=dt3;
145 Paso_Copy(n,u, u2);
146 }
147 }
148 if (errorCode == SOLVER_NO_ERROR) {
149 } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) {
150 /* if num_failures_max failures in a row: give up */
151 if (num_failures >= num_failures_max) {
152 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: No convergence after time step reductions.");
153 } else {
154 options->time_step_backtracking_used=TRUE;
155 if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n");
156 dt2=dt3*reduction_after_divergence_factor;
157 num_failures++;
158 Paso_Copy(n,u, u_save); /* reset initial value */
159 }
160 } else if (errorCode == SOLVER_INPUT_ERROR) {
161 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver.");
162 } else if (errorCode == SOLVER_MEMORY_ERROR) {
163 Esys_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allocation failed.");
164 } else if (errorCode == SOLVER_BREAKDOWN) {
165 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down.");
166 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
167 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm.");
168 } else {
169 Esys_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error.");
170 }
171 }
172 }
173 } /* end of time loop */
174 /*
175 * clean-up:
176 *
177 */
178
179 Paso_FCT_Solver_free(fctsolver);
180 Paso_ReactiveSolver_free(rsolver);
181 TMPMEMFREE(u_save);
182 TMPMEMFREE(u2);
183 }
184
185
186 double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp)
187 {
188 double dt_max=0., dt_R=LARGE_POSITIVE_FLOAT, dt_T=LARGE_POSITIVE_FLOAT;
189 index_t fail_loc=0;
190 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
191
192 if ( ! fctp->valid_matrices) {
193 /* set row-sum of mass_matrix */
194 Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);
195 {
196 /* check for positive entries in lumped_mass_matrix and set negative value for constraints */
197 index_t fail=0, i;
198 #pragma omp parallel private(i,fail_loc )
199 {
200
201 fail_loc=0;
202 #pragma omp for schedule(static)
203 for (i=0;i<n;++i) {
204 const double m_i=fctp->lumped_mass_matrix[i];
205 if ( m_i > 0 ) {
206 if (fctp->constraint_mask[i] > 0 ) fctp->lumped_mass_matrix[i]=-1.;
207 } else {
208 fail_loc=1;
209 }
210 }
211 #pragma omp critical
212 {
213 fail=MAX(fail, fail_loc);
214 }
215 }
216 #ifdef ESYS_MPI
217 {
218 fail_loc=fail;
219 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, fctp->mpi_info->comm);
220 }
221 #endif
222 if (fail < 0 )
223 Esys_setError(VALUE_ERROR, "Paso_FCTSolver_getSafeTimeStepSize: negative mass matrix entries detected.");
224 }
225 /* split off row-sum from transport_matrix */
226 Paso_SystemMatrix_makeZeroRowSums(fctp->transport_matrix,fctp->reactive_matrix);
227 /* get a copy of the main diagonal of the mass matrix */
228 Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix);
229
230 if (Esys_noError()) {
231 dt_R=Paso_ReactiveSolver_getSafeTimeStepSize(fctp);
232 dt_T=Paso_FCT_Solver_getSafeTimeStepSize(fctp);
233 }
234 if (Esys_noError()) {
235 fctp->dt_max_R=dt_R;
236 fctp->dt_max_T=dt_T;
237 fctp->valid_matrices=TRUE;
238 dt_max=MIN( 2 * dt_R,dt_T);
239 }
240 } else {
241 dt_max=MIN(2 * fctp->dt_max_R,fctp->dt_max_T); /* factor 2 as we use operator splitting */
242 }
243 return dt_max;
244 }

  ViewVC Help
Powered by ViewVC 1.1.26