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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4873 - (show annotations)
Wed Apr 16 06:38:51 2014 UTC (5 years, 5 months ago) by caltinay
File size: 9764 byte(s)
whitespace only changes.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 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 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************/
19
20 /* Paso: Transport solver
21 *
22 * solves Mu_t=Ku+q
23 *
24 * using an operator splitting approach (K=L+D where L is row sum zero and D
25 * is a diagonal matrix):
26 *
27 * - Mu_t=Du+q u(0)=u (reactive part)
28 * - Mv_t=Lv v(0)=u(dt/2) (transport part uses flux correction (FCT))
29 * - Mu_t=Du+q u(dt/2)=v(dt/2) (reactive part)
30 *
31 * to return u(dt)
32 *
33 */
34 /****************************************************************************/
35
36 /* Author: l.gross@uq.edu.au */
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 namespace paso {
47
48 void TransportProblem::solve(double* u, double dt, double* u0, double* q,
49 Options* options)
50 {
51 const double reduction_after_divergence_factor = 0.5;
52 const dim_t num_failures_max=50;
53
54 Performance pp;
55 ReactiveSolver* rsolver=NULL;
56 FCT_Solver* fctsolver=NULL;
57
58 dim_t i_substeps=0, n_substeps=1, num_failures=0;
59 double *u_save=NULL, *u2=NULL;
60 double dt2,t=0, dt3;
61 err_t errorCode=SOLVER_NO_ERROR;
62 const dim_t n = transport_matrix->getTotalNumRows();
63 options->time_step_backtracking_used = false;
64 options->num_iter=0;
65
66 if (dt <= 0.) {
67 Esys_setError(VALUE_ERROR, "TransportProblem::solve: dt must be positive.");
68 } else if (getBlockSize() > 1) {
69 Esys_setError(VALUE_ERROR, "TransportProblem::solve: block size >1 "
70 "is not supported.");
71 }
72 if (options->verbose) {
73 if (options->ode_solver == PASO_BACKWARD_EULER) {
74 printf("TransportProblem::solve: Backward Euler is used (dt = %e)\n",dt);
75 } else if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) {
76 printf("TransportProblem::solve: linear Crank-Nicolson is used (dt = %e).\n",dt);
77 } else if (options->ode_solver == PASO_CRANK_NICOLSON) {
78 printf("TransportProblem::solve: Crank-Nicolson is used (dt = %e).\n",dt);
79 } else {
80 Esys_setError(VALUE_ERROR, "TransportProblem::solve: unknown ODE solver.");
81 }
82 }
83 if (Esys_noError()) {
84 getSafeTimeStepSize();
85 // allocate memory
86 fctsolver = new FCT_Solver(shared_from_this(), options);
87 rsolver = new ReactiveSolver(shared_from_this());
88 u_save = new double[n];
89 u2 = new double[n];
90 }
91 if (Esys_noError()) {
92 // let the show begin!!!!
93 const double dt_R = dt_max_R;
94 const double dt_T = dt_max_T;
95 dt2 = dt;
96 if (dt_R < LARGE_POSITIVE_FLOAT)
97 dt2 = MIN(dt_R*2, dt); // as we half the step size for the RT bit
98 if (dt_T < LARGE_POSITIVE_FLOAT) {
99 if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON || options->ode_solver == PASO_CRANK_NICOLSON) {
100 dt2 = MIN(dt_T, dt);
101 } // PASO_BACKWARD_EULER does not require a restriction
102 }
103
104 num_failures = 0;
105 util::copy(n, u, u0); // copy initial value to return
106
107 while((dt-t) > dt*sqrt(EPSILON) && Esys_noError()) {
108 n_substeps = ceil((dt-t)/dt2);
109 if (n_substeps <= 0) {
110 Esys_setError(VALUE_ERROR, "TransportProblem::solve: time stepping break down.");
111 } else {
112 dt3 = (dt-t)/n_substeps;
113 if (options->verbose)
114 printf("TransportProblem::solve: number of substeps = %d "
115 "with dt = %e.\n", n_substeps, dt3);
116 // initialize the iteration matrix
117 fctsolver->initialize(dt3, options, &pp);
118 rsolver->initialize(dt3/2, options);
119 errorCode = SOLVER_NO_ERROR;
120
121 // start iteration
122 for (i_substeps=0; i_substeps<n_substeps &&
123 errorCode==SOLVER_NO_ERROR &&
124 Esys_noError(); i_substeps++) {
125 if (options->verbose)
126 printf("TransportProblem::solve: substep %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3);
127
128 // create copy for restart in case of failure
129 util::copy(n, u_save, u);
130 // update u
131
132 // Mu_t=Du+q u(0)=u
133 errorCode = rsolver->solve(u2, u, q, options, &pp);
134
135 // Mv_t=Lv v(0)=u(dt/2)
136 if (errorCode == SOLVER_NO_ERROR) {
137 errorCode = fctsolver->update(u, u2, options, &pp);
138
139 }
140 // Mu_t=Du+q u(dt/2)=v(dt/2)
141 if (errorCode == SOLVER_NO_ERROR) {
142 errorCode = rsolver->solve(u2, u, q, options, &pp);
143 }
144
145 if (errorCode == SOLVER_NO_ERROR) {
146 num_failures = 0;
147 t += dt3;
148 util::copy(n, u, u2);
149 }
150 }
151 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, "TransportProblem::solve: "
155 "No convergence after time step reductions.");
156 } else {
157 options->time_step_backtracking_used = true;
158 if (options->verbose)
159 printf("TransportProblem::solve: no convergence. "
160 "Time step size is reduced.\n");
161 dt2 = dt3*reduction_after_divergence_factor;
162 num_failures++;
163 util::copy(n, u, u_save); // reset initial value
164 }
165 } else if (errorCode == SOLVER_INPUT_ERROR) {
166 Esys_setError(VALUE_ERROR, "TransportProblem::solve: "
167 "input error for solver.");
168 } else if (errorCode == SOLVER_MEMORY_ERROR) {
169 Esys_setError(MEMORY_ERROR, "TransportProblem::solve: "
170 "memory allocation failed.");
171 } else if (errorCode == SOLVER_BREAKDOWN) {
172 Esys_setError(VALUE_ERROR, "TransportProblem::solve: "
173 "solver break down.");
174 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
175 Esys_setError(VALUE_ERROR, "TransportProblem::solve: "
176 "negative norm.");
177 } else if (errorCode != SOLVER_NO_ERROR) {
178 Esys_setError(SYSTEM_ERROR, "TransportProblem::solve: "
179 "general error.");
180 }
181 }
182 }
183 } // end of time loop
184
185 delete fctsolver;
186 delete rsolver;
187 delete[] u_save;
188 delete[] u2;
189 }
190
191 double TransportProblem::getSafeTimeStepSize()
192 {
193 double dt_max=0.;
194 const dim_t n = transport_matrix->getTotalNumRows();
195
196 if (!valid_matrices) {
197 // set row-sum of mass_matrix
198 mass_matrix->rowSum(lumped_mass_matrix);
199 // check for positive entries in lumped_mass_matrix and set
200 // negative value for constraints
201 int fail = 0;
202 #pragma omp parallel
203 {
204 index_t fail_loc = 0;
205 #pragma omp for
206 for (index_t i=0; i<n; ++i) {
207 const double m_i = lumped_mass_matrix[i];
208 if (m_i > 0) {
209 if (constraint_mask[i] > 0)
210 lumped_mass_matrix[i]=-1.;
211 } else {
212 fail_loc = 1;
213 }
214 }
215 #pragma omp critical
216 {
217 fail = MAX(fail, fail_loc);
218 }
219 }
220 #ifdef ESYS_MPI
221 int fail_loc = fail;
222 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, mpi_info->comm);
223 #endif
224 if (fail > 0)
225 Esys_setError(VALUE_ERROR, "TransportProblem::getSafeTimeStepSize: "
226 "negative mass matrix entries detected.");
227 // split off row-sum from transport_matrix
228 transport_matrix->makeZeroRowSums(reactive_matrix);
229 // get a copy of the main diagonal of the mass matrix
230 mass_matrix->copyFromMainDiagonal(main_diagonal_mass_matrix);
231
232 if (Esys_noError()) {
233 const double dt_R = ReactiveSolver::getSafeTimeStepSize(shared_from_this());
234 const double dt_T = FCT_Solver::getSafeTimeStepSize(shared_from_this());
235 dt_max_R = dt_R;
236 dt_max_T = dt_T;
237 valid_matrices = true;
238 dt_max = MIN(2*dt_R, dt_T);
239 }
240 } else {
241 // factor 2 as we use operator splitting
242 dt_max = MIN(2*dt_max_R, dt_max_T);
243 }
244 return dt_max;
245 }
246
247 } // namespace paso
248

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/Transport_solve.cpp:3531-3826 /branches/lapack2681/paso/src/Transport_solve.cpp:2682-2741 /branches/pasowrap/paso/src/Transport_solve.cpp:3661-3674 /branches/py3_attempt2/paso/src/Transport_solve.cpp:3871-3891 /branches/restext/paso/src/Transport_solve.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/Transport_solve.cpp:3669-3791 /branches/stage3.0/paso/src/Transport_solve.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/Transport_solve.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/Transport_solve.cpp:3517-3974 /release/3.0/paso/src/Transport_solve.cpp:2591-2601 /trunk/paso/src/Transport_solve.cpp:4257-4344 /trunk/ripley/test/python/paso/src/Transport_solve.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26