/[escript]/branches/doubleplusgood/paso/src/Transport_solve.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/Transport_solve.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4324 - (show annotations)
Wed Mar 20 00:55:44 2013 UTC (5 years, 10 months ago) by jfenwick
File size: 9314 byte(s)
continues
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Paso: Transport solver
20 *
21 * solves Mu_t=Ku+q
22 *
23 * using an operator splitting approach (K=L+D where L is row sum zero and D is a diagonal matrix):
24 *
25 * - Mu_t=Du+q u(0)=u (reactive part)
26 * - Mv_t=Lv v(0)=u(dt/2) (transport part uses flux correction (FCT))
27 * - Mu_t=Du+q u(dt/2)=v(dt/2) (reactive part)
28 *
29 * to return u(dt)
30 *
31 */
32 /************************************************************************************/
33
34 /* Author: l.gross@uq.edu.au */
35
36 /************************************************************************************/
37
38
39
40 #include "Transport.h"
41 #include "FCT_Solver.h"
42 #include "ReactiveSolver.h"
43 #include "Solver.h"
44 #include "PasoUtil.h"
45
46
47 void Paso_TransportProblem_solve(Paso_TransportProblem* fctp, double* u, double dt, double* u0, double* q, Paso_Options* options) {
48 const double reduction_after_divergence_factor = 0.5;
49 const dim_t num_failures_max=50;
50
51 Paso_Performance pp;
52 Paso_ReactiveSolver *rsolver=NULL;
53 Paso_FCT_Solver *fctsolver=NULL;
54
55
56 dim_t i_substeps=0, n_substeps=1, num_failures=0;
57 double *u_save=NULL, *u2=NULL;
58 double dt2,t=0, dt3;
59 err_t errorCode=SOLVER_NO_ERROR;
60 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
61 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);
62 options->time_step_backtracking_used=FALSE;
63 options->num_iter=0;
64
65 if (dt<=0.) {
66 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");
67 }
68 if (blockSize>1) {
69 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >1 is not supported.");
70 }
71 if (options->verbose) {
72 if (options->ode_solver == PASO_BACKWARD_EULER) {
73 printf("Paso_TransportProblem_solve: Backward Euler is used (dt = %e)\n",dt);
74 } else if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) {
75 printf("Paso_TransportProblem_solve: linear Crank-Nicolson is used (dt = %e).\n",dt);
76 } else if (options->ode_solver == PASO_CRANK_NICOLSON) {
77 printf("Paso_TransportProblem_solve: Crank-Nicolson is used (dt = %e).\n",dt);
78 } else {
79 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: unknown ODE solver.");
80 }
81 }
82 if (Esys_noError()) {
83 Paso_TransportProblem_getSafeTimeStepSize(fctp);
84 /*
85 * allocate memory
86 *
87 */
88 fctsolver=Paso_FCT_Solver_alloc(fctp,options);
89 rsolver=Paso_ReactiveSolver_alloc(fctp);
90 u_save=new double[n];
91 u2=new double[n];
92 Esys_checkPtr(u_save);
93 Esys_checkPtr(u2);
94 }
95 if (Esys_noError()) {
96 /*
97 * let the show begin!!!!
98 *
99 */
100 const double dt_R=fctp->dt_max_R;
101 const double dt_T=fctp->dt_max_T;
102 dt2=dt;
103 if (dt_R < LARGE_POSITIVE_FLOAT) dt2 = MIN(dt_R *2, dt); /* as we half the steps size for the RT bit */
104 if (dt_T < LARGE_POSITIVE_FLOAT) {
105 if ( ( options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) || ( options->ode_solver == PASO_CRANK_NICOLSON ) ) {
106 dt2 = MIN(dt_T, dt);
107 } /* PASO_BACKWARD_EULER does not require a restriction */
108 }
109
110 num_failures=0;
111 Paso_Copy(n,u,u0); /* copy initial value to return */
112
113 while( (dt-t)>dt*sqrt(EPSILON) && Esys_noError()) {
114
115 n_substeps=ceil((dt-t)/dt2);
116 if ( n_substeps <= 0) {
117 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: time stepping break down.");
118 } else {
119 dt3=(dt-t)/n_substeps;
120 if (options->verbose) printf("Paso_TransportProblem_solve: number of substeps = %d with dt = %e.\n",n_substeps,dt3);
121 /* initialize the iteration matrix */
122 Paso_FCT_Solver_initialize(dt3, fctsolver, options, &pp);
123 Paso_ReactiveSolver_initialize(dt3/2, rsolver, options);
124 errorCode = SOLVER_NO_ERROR;
125
126 /* start iteration */
127 for (i_substeps =0; (i_substeps<n_substeps) && (errorCode==SOLVER_NO_ERROR) && Esys_noError(); i_substeps++) {
128 if (options->verbose) printf("Paso_TransportProblem_solve: substep %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3);
129 Paso_Copy(n,u_save, u); /* create copy for restart in case of failure */
130 /* updates u */
131 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u2, u ,q, options, &pp); /* Mu_t=Du+q u(0)=u */
132 if (errorCode == SOLVER_NO_ERROR) {
133 errorCode=Paso_FCT_Solver_update(fctsolver, u, u2,options, &pp); /* Mv_t=Lv v(0)=u(dt/2) */
134
135 }
136 if (errorCode == SOLVER_NO_ERROR) {
137 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u2, u, q, options, &pp); /* Mu_t=Du+q u(dt/2)=v(dt/2) */
138 }
139
140
141 /*
142 * error handling:
143 */
144 if (errorCode == SOLVER_NO_ERROR) {
145 num_failures=0;
146 t+=dt3;
147 Paso_Copy(n,u, u2);
148 }
149 }
150 if (errorCode == SOLVER_NO_ERROR) {
151 } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) {
152 /* if num_failures_max failures in a row: give up */
153 if (num_failures >= num_failures_max) {
154 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: No convergence after time step reductions.");
155 } else {
156 options->time_step_backtracking_used=TRUE;
157 if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n");
158 dt2=dt3*reduction_after_divergence_factor;
159 num_failures++;
160 Paso_Copy(n,u, u_save); /* reset initial value */
161 }
162 } else if (errorCode == SOLVER_INPUT_ERROR) {
163 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver.");
164 } else if (errorCode == SOLVER_MEMORY_ERROR) {
165 Esys_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allocation failed.");
166 } else if (errorCode == SOLVER_BREAKDOWN) {
167 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down.");
168 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
169 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm.");
170 } else {
171 Esys_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error.");
172 }
173 }
174 }
175 } /* end of time loop */
176 /*
177 * clean-up:
178 *
179 */
180
181 Paso_FCT_Solver_free(fctsolver);
182 Paso_ReactiveSolver_free(rsolver);
183 delete[] u_save;
184 delete[] u2;
185 }
186
187
188 double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp)
189 {
190 double dt_max=0., dt_R=LARGE_POSITIVE_FLOAT, dt_T=LARGE_POSITIVE_FLOAT;
191 index_t fail_loc=0;
192 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
193
194 if ( ! fctp->valid_matrices) {
195 /* set row-sum of mass_matrix */
196 Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);
197 {
198 /* check for positive entries in lumped_mass_matrix and set negative value for constraints */
199 index_t fail=0, i;
200 #pragma omp parallel private(i,fail_loc )
201 {
202
203 fail_loc=0;
204 #pragma omp for schedule(static)
205 for (i=0;i<n;++i) {
206 const double m_i=fctp->lumped_mass_matrix[i];
207 if ( m_i > 0 ) {
208 if (fctp->constraint_mask[i] > 0 ) fctp->lumped_mass_matrix[i]=-1.;
209 } else {
210 fail_loc=1;
211 }
212 }
213 #pragma omp critical
214 {
215 fail=MAX(fail, fail_loc);
216 }
217 }
218 #ifdef ESYS_MPI
219 {
220 fail_loc=fail;
221 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, fctp->mpi_info->comm);
222 }
223 #endif
224 if (fail < 0 )
225 Esys_setError(VALUE_ERROR, "Paso_FCTSolver_getSafeTimeStepSize: negative mass matrix entries detected.");
226 }
227 /* split off row-sum from transport_matrix */
228 Paso_SystemMatrix_makeZeroRowSums(fctp->transport_matrix,fctp->reactive_matrix);
229 /* get a copy of the main diagonal of the mass matrix */
230 Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix);
231
232 if (Esys_noError()) {
233 dt_R=Paso_ReactiveSolver_getSafeTimeStepSize(fctp);
234 dt_T=Paso_FCT_Solver_getSafeTimeStepSize(fctp);
235 }
236 if (Esys_noError()) {
237 fctp->dt_max_R=dt_R;
238 fctp->dt_max_T=dt_T;
239 fctp->valid_matrices=TRUE;
240 dt_max=MIN( 2 * dt_R,dt_T);
241 }
242 } else {
243 dt_max=MIN(2 * fctp->dt_max_R,fctp->dt_max_T); /* factor 2 as we use operator splitting */
244 }
245 return dt_max;
246 }

  ViewVC Help
Powered by ViewVC 1.1.26