/[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 6000 - (show annotations)
Tue Mar 1 00:24:43 2016 UTC (3 years, 1 month ago) by caltinay
File size: 9226 byte(s)
a few more include rearrangements.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The 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 #include <iostream>
47
48 namespace paso {
49
50 void TransportProblem::solve(double* u, double dt, double* u0, double* q,
51 Options* options)
52 {
53 const real_t EPSILON = escript::DataTypes::real_t_eps();
54 const real_t LARGE_POSITIVE_FLOAT = escript::DataTypes::real_t_max();
55 const double reduction_after_divergence_factor = 0.5;
56 const dim_t num_failures_max=50;
57
58 Performance pp;
59 ReactiveSolver* rsolver=NULL;
60 FCT_Solver* fctsolver=NULL;
61
62 dim_t i_substeps=0, n_substeps=1, num_failures=0;
63 double *u_save=NULL, *u2=NULL;
64 double dt2,t=0, dt3;
65 SolverResult errorCode=NoError;
66 const dim_t n = transport_matrix->getTotalNumRows();
67 options->time_step_backtracking_used = false;
68 options->num_iter=0;
69
70 if (dt <= 0.) {
71 throw PasoException("TransportProblem::solve: dt must be positive.");
72 } else if (getBlockSize() > 1) {
73 throw PasoException("TransportProblem::solve: block size >1 "
74 "is not supported.");
75 }
76 if (options->verbose) {
77 if (options->ode_solver == PASO_BACKWARD_EULER) {
78 printf("TransportProblem::solve: Backward Euler is used (dt = %e)\n",dt);
79 } else if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) {
80 printf("TransportProblem::solve: linear Crank-Nicolson is used (dt = %e).\n",dt);
81 } else if (options->ode_solver == PASO_CRANK_NICOLSON) {
82 printf("TransportProblem::solve: Crank-Nicolson is used (dt = %e).\n",dt);
83 } else {
84 throw PasoException("TransportProblem::solve: unknown ODE solver.");
85 }
86 }
87 getSafeTimeStepSize();
88 // allocate memory
89 fctsolver = new FCT_Solver(shared_from_this(), options);
90 rsolver = new ReactiveSolver(shared_from_this());
91 u_save = new double[n];
92 u2 = new double[n];
93
94 // let the show begin!!!!
95 const double dt_R = dt_max_R;
96 const double dt_T = dt_max_T;
97 dt2 = dt;
98 if (dt_R < LARGE_POSITIVE_FLOAT)
99 dt2 = std::min(dt_R*2, dt); // as we half the step size for the RT bit
100 if (dt_T < LARGE_POSITIVE_FLOAT) {
101 if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON || options->ode_solver == PASO_CRANK_NICOLSON) {
102 dt2 = std::min(dt_T, dt);
103 } // PASO_BACKWARD_EULER does not require a restriction
104 }
105
106 num_failures = 0;
107 util::copy(n, u, u0); // copy initial value to return
108
109 while((dt-t) > dt*sqrt(EPSILON)) {
110 n_substeps = ceil((dt-t)/dt2);
111 if (n_substeps <= 0) {
112 throw PasoException("TransportProblem::solve: time stepping break down.");
113 } else {
114 dt3 = (dt-t)/n_substeps;
115 if (options->verbose) {
116 std::cout << "TransportProblem::solve: number of substeps = "
117 << n_substeps << " with dt = " << dt3 << "."
118 << std::endl;
119 }
120 // initialize the iteration matrix
121 fctsolver->initialize(dt3, options, &pp);
122 rsolver->initialize(dt3/2, options);
123 errorCode = NoError;
124
125 // start iteration
126 for (i_substeps=0; i_substeps<n_substeps &&
127 errorCode==NoError; i_substeps++) {
128 if (options->verbose) {
129 std::cout << "TransportProblem::solve: substep "
130 << i_substeps << " of " << n_substeps << " at t = "
131 << (t+dt3) << " (dt = " << dt3 << ")" << std::endl;
132 }
133
134 // create copy for restart in case of failure
135 util::copy(n, u_save, u);
136 // update u
137
138 // Mu_t=Du+q u(0)=u
139 errorCode = rsolver->solve(u2, u, q, options, &pp);
140
141 // Mv_t=Lv v(0)=u(dt/2)
142 if (errorCode == NoError) {
143 errorCode = fctsolver->update(u, u2, options, &pp);
144
145 }
146 // Mu_t=Du+q u(dt/2)=v(dt/2)
147 if (errorCode == NoError) {
148 errorCode = rsolver->solve(u2, u, q, options, &pp);
149 }
150
151 if (errorCode == NoError) {
152 num_failures = 0;
153 t += dt3;
154 util::copy(n, u, u2);
155 }
156 }
157 if (errorCode == MaxIterReached || errorCode == Divergence) {
158 // if num_failures_max failures in a row: give up
159 if (num_failures >= num_failures_max) {
160 throw PasoException("TransportProblem::solve: "
161 "No convergence after time step reductions.");
162 } else {
163 options->time_step_backtracking_used = true;
164 if (options->verbose) {
165 std::cout << "TransportProblem::solve: "
166 << "no convergence. Time step size is reduced."
167 << std::endl;
168 }
169 dt2 = dt3*reduction_after_divergence_factor;
170 num_failures++;
171 util::copy(n, u, u_save); // reset initial value
172 }
173 } else if (errorCode == InputError) {
174 throw PasoException("TransportProblem::solve: input error for solver.");
175 } else if (errorCode == MemoryError) {
176 throw PasoException("TransportProblem::solve: memory allocation failed.");
177 } else if (errorCode == Breakdown) {
178 throw PasoException("TransportProblem::solve: solver break down.");
179 } else if (errorCode == NegativeNormError) {
180 throw PasoException("TransportProblem::solve: negative norm.");
181 } else if (errorCode != NoError) {
182 throw PasoException("TransportProblem::solve: general error.");
183 }
184 }
185 } // end of time loop
186
187 delete fctsolver;
188 delete rsolver;
189 delete[] u_save;
190 delete[] u2;
191 }
192
193 double TransportProblem::getSafeTimeStepSize() const
194 {
195 double dt_max=0.;
196 const dim_t n = transport_matrix->getTotalNumRows();
197
198 if (!valid_matrices) {
199 // set row-sum of mass_matrix
200 mass_matrix->rowSum(lumped_mass_matrix);
201 // check for positive entries in lumped_mass_matrix and set
202 // negative value for constraints
203 int fail = 0;
204 #pragma omp parallel
205 {
206 index_t fail_loc = 0;
207 #pragma omp for
208 for (index_t i=0; i<n; ++i) {
209 const double m_i = lumped_mass_matrix[i];
210 if (m_i > 0) {
211 if (constraint_mask[i] > 0)
212 lumped_mass_matrix[i]=-1.;
213 } else {
214 fail_loc = 1;
215 }
216 }
217 #pragma omp critical
218 {
219 fail = std::max(fail, fail_loc);
220 }
221 }
222 #ifdef ESYS_MPI
223 int fail_loc = fail;
224 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, mpi_info->comm);
225 #endif
226 if (fail > 0)
227 throw PasoException("TransportProblem::getSafeTimeStepSize: "
228 "negative mass matrix entries detected.");
229 // split off row-sum from transport_matrix
230 transport_matrix->makeZeroRowSums(reactive_matrix);
231 // get a copy of the main diagonal of the mass matrix
232 mass_matrix->copyFromMainDiagonal(main_diagonal_mass_matrix);
233
234 const double dt_R = ReactiveSolver::getSafeTimeStepSize(shared_from_this());
235 const double dt_T = FCT_Solver::getSafeTimeStepSize(shared_from_this());
236 dt_max_R = dt_R;
237 dt_max_T = dt_T;
238 valid_matrices = true;
239 dt_max = std::min(2*dt_R, dt_T);
240 } else {
241 // factor 2 as we use operator splitting
242 dt_max = std::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/4.0fordebian/paso/src/Transport_solve.cpp:5567-5588 /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 /release/4.0/paso/src/Transport_solve.cpp:5380-5406 /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