112 template <
class Grid ,
116 void solve(
const Grid& g ,
117 const SourceTerms* src ,
121 LinearSolver& linsolve,
124 typedef typename JacobianSystem::vector_type vector_type;
125 typedef typename JacobianSystem::matrix_type matrix_type;
127 asm_.createSystem(g, sys_);
128 model_.initStep(state, g, sys_);
129 model_.initIteration(state, g, sys_);
131 MZero<matrix_type>::zero(sys_.writableMatrix());
132 VZero<vector_type>::zero(sys_.vector().writableResidual());
134 asm_.assemble(state, g, src, dt, sys_);
136 const double nrm_res0 =
137 VNorm<vector_type>::norm(sys_.vector().residual());
139 rpt.norm_res = nrm_res0;
146 VZero<vector_type>::zero(sys_.vector().writableIncrement());
148 linsolve.solve(sys_.matrix(),
149 sys_.vector().residual(),
150 sys_.vector().writableIncrement());
152 VNeg<vector_type>::negate(sys_.vector().writableIncrement());
157 VNorm<vector_type>::norm(sys_.vector().increment());
162 bool finished=rpt.norm_res<ctrl.atol;
165 vector_type dx_old(sys_.vector().increment());
166 vector_type x_old(sys_.vector().solution());
169 VAsgn<vector_type>::assign(alpha, dx_old,
170 sys_.vector().writableIncrement());
171 VAsgn<vector_type>::assign(x_old,
172 sys_.vector().writableSolution());
174 sys_.vector().addIncrement();
175 bool init = model_.initIteration(state, g, sys_);
177 MZero<matrix_type>::zero(sys_.writableMatrix());
178 VZero<vector_type>::zero(sys_.vector().writableResidual());
179 asm_.assemble(state, g, src, dt, sys_);
180 residual = VNorm<vector_type>::norm(sys_.vector().residual());
181 if (ctrl.verbosity > 1){
182 std::cout <<
"Line search iteration " << std::scientific << lin_it
183 <<
" norm :" << residual <<
" alpha " << alpha <<
'\n';
186 if (ctrl.verbosity > 1){
187 std::cout <<
"Line search iteration " << std::scientific << lin_it
188 <<
" Value out of range, continue search. alpha " << alpha <<
'\n';
192 finished=(residual < rpt.norm_res) || (lin_it> ctrl.max_it_ls);
196 VNorm<vector_type>::norm(sys_.vector().residual());
199 if (ctrl.verbosity > 0){
200 std::cout <<
"Iteration " << std::scientific << rpt.nit
201 <<
" norm :" << rpt.norm_res <<
" alpha " << alpha << std::endl;
203 done = (rpt.norm_res < ctrl.atol) ||
204 (rpt.norm_res < ctrl.rtol * nrm_res0) ||
205 (rpt.norm_dx < ctrl.dxtol) ||
206 (lin_it > ctrl.max_it_ls) ||
207 (rpt.nit == ctrl.max_it);
210 model_.finishStep(g, sys_.vector().solution(), state);
212 if (rpt.norm_res < ctrl.atol) { rpt.flag = 1; }
213 else if (rpt.norm_res < ctrl.rtol * nrm_res0) { rpt.flag = 2; }
214 else if (rpt.norm_dx < ctrl.dxtol) { rpt.flag = 3; }
215 else { rpt.flag = -1; }
223 using Model::initStep;
224 using Model::initIteration;
225 using Model::finishIteration;
226 using Model::finishStep;