← Back to branch summary

~ubuntu-branches/ubuntu/wily/dolfin/wily-proposed

« back to all changes in this revision

Viewing changes to dolfin/multistage/PointIntegralSolver.cpp

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2014-09-22 14:35:34 UTC
  • mfrom: (1.1.17) (19.1.23 sid)
  • Revision ID: package-import@ubuntu.com-20140922143534-0yi89jyuqbgdxwm9
Tags: 1.4.0+dfsg-4
http://bugs.debian.org/760658
* debian/control: Disable libcgal-dev on i386, mipsel and sparc.
* debian/rules: Remove bad directives in pkg-config file dolfin.pc
  (closes: #760658).
* Remove debian/libdolfin-dev.lintian-overrides.

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/multistage/PointIntegralSolver.cpp
16
16
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
17
17
//
18
18
// First added:  2013-02-15
19
 
// Last changed: 2013-05-13
 
19
// Last changed: 2014-04-09
20
20
 
21
21
#include <cmath>
22
 
#include <boost/make_shared.hpp>
23
 
#include <Eigen/Dense>
 
22
#include <algorithm>
 
23
#include <memory>
24
24
 
25
25
#include <dolfin/log/log.h>
26
26
#include <dolfin/common/Timer.h>
43
43
using namespace dolfin;
44
44
 
45
45
//-----------------------------------------------------------------------------
46
 
PointIntegralSolver::PointIntegralSolver(boost::shared_ptr<MultiStageScheme> scheme) :
47
 
  Variable("PointIntegralSolver", "unamed"),
48
 
  _scheme(scheme), _vertex_map(), _ufcs(), _coefficient_index(),
49
 
  _retabulate_J(true)
 
46
PointIntegralSolver::PointIntegralSolver(std::shared_ptr<MultiStageScheme> scheme) :
 
47
  Variable("PointIntegralSolver", "unnamed"), _scheme(scheme),
 
48
  _mesh(_scheme->stage_forms()[0][0]->mesh()),
 
49
  _dofmap(*_scheme->stage_forms()[0][0]->function_space(0)->dofmap()),
 
50
  _system_size(_dofmap.num_entity_dofs(0)), _dof_offset(_mesh.type().num_entities(0)),
 
51
  _num_stages(_scheme->stage_forms().size()), _local_to_local_dofs(_system_size),
 
52
  _vertex_map(), _local_to_global_dofs(_system_size),
 
53
  _local_stage_solutions(_scheme->stage_solutions().size()),
 
54
  _u0(_system_size), _residual(_system_size), _y(_system_size), _dx(_system_size),
 
55
  _ufcs(), _coefficient_index(), _recompute_jacobian(),
 
56
  _jacobians(), _eta(1.0), _num_jacobian_computations(0)
50
57
{
 
58
  Timer construct_pis("Construct PointIntegralSolver");
 
59
 
51
60
  // Set parameters
52
61
  parameters = default_parameters();
53
62
 
55
64
  _init();
56
65
}
57
66
//-----------------------------------------------------------------------------
 
67
PointIntegralSolver::~PointIntegralSolver()
 
68
{
 
69
  // Do nothing
 
70
}
 
71
//-----------------------------------------------------------------------------
 
72
void PointIntegralSolver::reset_newton_solver()
 
73
{
 
74
  _eta = parameters("newton_solver")["eta_0"];
 
75
 
 
76
  for (unsigned int i=0; i < _recompute_jacobian.size(); i++)
 
77
    _recompute_jacobian[i] = true;
 
78
}
 
79
//-----------------------------------------------------------------------------
 
80
void PointIntegralSolver::reset_stage_solutions()
 
81
{
 
82
  // Iterate over stage forms
 
83
  for (unsigned int stage=0; stage<_num_stages; stage++)
 
84
  {
 
85
    // Reset global stage solutions
 
86
    *_scheme->stage_solutions()[stage]->vector() = 0.0;
 
87
 
 
88
    // Reset local stage solutions
 
89
    for (unsigned int row=0; row < _system_size; row++)
 
90
      _local_stage_solutions[stage][row] = 0.0;
 
91
 
 
92
  }
 
93
}
 
94
//-----------------------------------------------------------------------------
58
95
void PointIntegralSolver::step(double dt)
59
96
{
 
97
 
 
98
  const bool reset_stage_solutions_ = parameters["reset_stage_solutions"];
 
99
  const bool reset_newton_solver_ = parameters("newton_solver")["reset_each_step"];
 
100
 
 
101
  // Check for reseting stage solutions
 
102
  if (reset_stage_solutions_)
 
103
    reset_stage_solutions();
 
104
 
 
105
  // Check for reseting newtonsolver for each time step
 
106
  if (reset_newton_solver_)
 
107
    reset_newton_solver();
 
108
 
 
109
  Timer t_step("PointIntegralSolver::step");
 
110
 
60
111
  dolfin_assert(dt > 0.0);
61
112
 
62
 
  Timer t_step_stage("Step: set the stage");
63
 
 
64
113
  // Update time constant of scheme
65
114
  *_scheme->dt() = dt;
66
115
 
67
116
  // Time at start of timestep
68
117
  const double t0 = *_scheme->t();
69
118
 
70
 
  // Extract mesh
71
 
  const Mesh& mesh = _scheme->stage_forms()[0][0]->mesh();
72
 
 
73
 
  // Collect ref to dof map only need one as we require same trial and
74
 
  // test space for all forms
75
 
  const GenericDofMap& dofmap
76
 
    = *_scheme->stage_forms()[0][0]->function_space(0)->dofmap();
77
 
 
78
 
  // Get size of system (num dofs per vertex)
79
 
  const unsigned int N = dofmap.num_entity_dofs(0);
80
 
  const unsigned int dof_offset = mesh.type().num_entities(0);
81
 
  const unsigned int num_stages = _scheme->stage_forms().size();
82
 
 
83
 
  // Local solution vector at start of time step
84
 
 
85
 
  Eigen::VectorXd u0(N);
86
 
 
87
 
  // Local stage solutions
88
 
  std::vector<Eigen::VectorXd >
89
 
    local_stage_solutions(_scheme->stage_solutions().size());
90
 
  for (unsigned int stage = 0; stage < num_stages; stage++)
91
 
    local_stage_solutions[stage].resize(N);
92
 
 
93
 
  // Update off-process coefficients
94
 
  for (unsigned int i = 0; i < num_stages; i++)
95
 
  {
96
 
    for (unsigned int j = 0; j < _scheme->stage_forms()[i].size(); j++)
97
 
    {
98
 
      const std::vector<boost::shared_ptr<const GenericFunction> >
99
 
        coefficients = _scheme->stage_forms()[i][j]->coefficients();
100
 
 
101
 
      for (unsigned int k = 0; k < coefficients.size(); ++k)
102
 
        coefficients[k]->update();
103
 
    }
104
 
  }
105
 
 
106
 
  // Local to local dofs to be used in tabulate entity dofs
107
 
  std::vector<std::size_t> local_to_local_dofs(N);
108
 
 
109
 
  // Local to global dofs used when solution is fanned out to global vector
110
 
  std::vector<dolfin::la_index> local_to_global_dofs(N);
 
119
  // Get ownership range
 
120
  const dolfin::la_index n0 = _dofmap.ownership_range().first;
 
121
  const dolfin::la_index n1 = _dofmap.ownership_range().second;
111
122
 
112
123
  // Iterate over vertices
113
 
  //Progress p("Solving local point integral problems", mesh.num_vertices());
114
 
 
115
 
  t_step_stage.stop();
116
 
 
117
124
  ufc::cell ufc_cell;
118
125
  std::vector<double> vertex_coordinates;
119
 
  for (std::size_t vert_ind = 0; vert_ind < mesh.num_vertices(); ++vert_ind)
 
126
  for (std::size_t vert_ind=0; vert_ind < _mesh.num_vertices(); ++vert_ind)
120
127
  {
121
 
    Timer t_vert("Step: update vert");
122
128
 
123
129
    // Cell containing vertex
124
 
    const Cell cell(mesh, _vertex_map[vert_ind].first);
 
130
    const Cell cell(_mesh, _vertex_map[vert_ind].first);
125
131
    cell.get_vertex_coordinates(vertex_coordinates);
126
132
    cell.get_cell_data(ufc_cell);
127
133
 
128
134
    // Get all dofs for cell
129
135
    // FIXME: Shold we include logics about empty dofmaps?
130
136
    const std::vector<dolfin::la_index>& cell_dofs
131
 
      = dofmap.cell_dofs(cell.index());
132
 
 
133
 
    // Local vertex ind
134
 
    const unsigned int local_vert = _vertex_map[vert_ind].second;
 
137
      = _dofmap.cell_dofs(cell.index());
135
138
 
136
139
    // Tabulate local-local dofmap
137
 
    dofmap.tabulate_entity_dofs(local_to_local_dofs, 0, local_vert);
138
 
 
139
 
    // Fill local to global dof map
140
 
    for (unsigned int row = 0; row < N; row++)
141
 
      local_to_global_dofs[row] = cell_dofs[local_to_local_dofs[row]];
142
 
 
143
 
    t_vert.stop();
 
140
    _dofmap.tabulate_entity_dofs(_local_to_local_dofs, 0,
 
141
                                 _vertex_map[vert_ind].second);
 
142
 
 
143
    // Fill local to global dof map and check that the dof is owned
 
144
    bool owns_all_dofs = true;
 
145
    for (unsigned int row=0; row < _system_size; row++)
 
146
    {
 
147
      _local_to_global_dofs[row] = cell_dofs[_local_to_local_dofs[row]];
 
148
      if (_local_to_global_dofs[row] < n0 || n1 <=_local_to_global_dofs[row])
 
149
      {
 
150
        owns_all_dofs = false;
 
151
        break;
 
152
      }
 
153
    }
 
154
 
 
155
    // If not owning all dofs
 
156
    if (!owns_all_dofs)
 
157
      continue;
144
158
 
145
159
    // Iterate over stage forms
146
 
    for (unsigned int stage = 0; stage < num_stages; stage++)
 
160
    for (unsigned int stage=0; stage<_num_stages; stage++)
147
161
    {
148
 
      // Update time
149
 
      *_scheme->t() = t0 + dt*_scheme->dt_stage_offset()[stage];
 
162
 
 
163
      // Update cell
 
164
      // TODO: Pass suitable bool vector here to avoid tabulating all coefficient dofs:
 
165
      _ufcs[stage][0]->update(cell, vertex_coordinates, ufc_cell);
 
166
                              //some_integral.enabled_coefficients());
150
167
 
151
168
      // Check if we have an explicit stage (only 1 form)
152
169
      if (_ufcs[stage].size() == 1)
153
170
      {
154
 
        Timer t_expl("Explicit stage");
155
 
 
156
 
        // Point integral
157
 
        const ufc::point_integral& integral
158
 
          = *_ufcs[stage][0]->default_point_integral;
159
 
 
160
 
        // Update to current cell
161
 
        Timer t_expl_update("Explicit stage: update_cell");
162
 
        _ufcs[stage][0]->update(cell, vertex_coordinates, ufc_cell);
163
 
        t_expl_update.stop();
164
 
 
165
 
        // Tabulate cell tensor
166
 
        Timer t_expl_tt("Explicit stage: tabulate_tensor");
167
 
        integral.tabulate_tensor(_ufcs[stage][0]->A.data(),
168
 
                                 _ufcs[stage][0]->w(),
169
 
                                 vertex_coordinates.data(),
170
 
                                 local_vert);
171
 
        t_expl_tt.stop();
172
 
 
173
 
        // Extract vertex dofs from tabulated tensor and put them into the local
174
 
        // stage solution vector
175
 
        // Extract vertex dofs from tabulated tensor
176
 
        for (unsigned int row=0; row < N; row++)
177
 
        {
178
 
          local_stage_solutions[stage](row)
179
 
            = _ufcs[stage][0]->A[local_to_local_dofs[row]];
180
 
        }
181
 
 
182
 
        // Put solution back into global stage solution vector
183
 
        Timer t_expl_set("Explicit stage: set");
184
 
        _scheme->stage_solutions()[stage]->vector()->set(
185
 
          local_stage_solutions[stage].data(), N,
186
 
          local_to_global_dofs.data());
 
171
        _solve_explicit_stage(vert_ind, stage, ufc_cell, vertex_coordinates);
187
172
      }
188
173
 
189
174
      // or an implicit stage (2 forms)
190
175
      else
191
176
      {
192
 
        Timer t_impl("Implicit stage");
193
 
 
194
 
        const Parameters& newton_parameters = parameters("newton_solver");
195
 
 
196
 
        // Local solution
197
 
        Eigen::VectorXd& u = local_stage_solutions[stage];
198
 
 
199
 
        unsigned int newton_iteration = 0;
200
 
        bool newton_converged = false;
201
 
        bool jacobian_retabulated = false;
202
 
        const std::size_t maxiter = newton_parameters["maximum_iterations"];
203
 
        const bool reuse_jacobian = newton_parameters["reuse_jacobian"];
204
 
        const std::size_t iterations_to_retabulate_jacobian
205
 
          = newton_parameters["iterations_to_retabulate_jacobian"];
206
 
        const double relaxation = newton_parameters["relaxation_parameter"];
207
 
        const std::string convergence_criterion
208
 
          = newton_parameters["convergence_criterion"];
209
 
        const double rtol = newton_parameters["relative_tolerance"];
210
 
        const double atol = newton_parameters["absolute_tolerance"];
211
 
        const bool report = newton_parameters["report"];
212
 
        //const int num_threads = 0;
213
 
 
214
 
        /// Most recent residual and intitial residual
215
 
        double residual = 1.0;
216
 
        double residual0 = 1.0;
217
 
        double relative_residual = 1.0;
218
 
 
219
 
        //const double relaxation = 1.0;
220
 
 
221
 
        // Initialize la structures
222
 
        Eigen::VectorXd F(N), dx(N);
223
 
 
224
 
        // Get point integrals
225
 
        const ufc::point_integral& F_integral
226
 
          = *_ufcs[stage][0]->default_point_integral;
227
 
        const ufc::point_integral& J_integral
228
 
          = *_ufcs[stage][1]->default_point_integral;
229
 
 
230
 
        // Update to current cell. This only need to be done once for
231
 
        // each stage and vertex
232
 
        Timer t_impl_update("Implicit stage: update_cell");
233
 
        _ufcs[stage][0]->update(cell, vertex_coordinates, ufc_cell);
234
 
        _ufcs[stage][1]->update(cell, vertex_coordinates, ufc_cell);
235
 
        t_impl_update.stop();
236
 
 
237
 
        // Tabulate an initial residual solution
238
 
        Timer t_impl_tt_F("Implicit stage: tabulate_tensor (F)");
239
 
        F_integral.tabulate_tensor(_ufcs[stage][0]->A.data(),
240
 
                                   _ufcs[stage][0]->w(),
241
 
                                   vertex_coordinates.data(),
242
 
                                   local_vert);
243
 
        t_impl_tt_F.stop();
244
 
 
245
 
        // Extract vertex dofs from tabulated tensor, together with
246
 
        // the old stage solution
247
 
        Timer t_impl_update_F("Implicit stage: update_F");
248
 
        for (unsigned int row = 0; row < N; row++)
249
 
        {
250
 
          F(row) = _ufcs[stage][0]->A[local_to_local_dofs[row]];
251
 
 
252
 
          // Grab old value of stage solution as an initial start
253
 
          // value. This value was also used to tabulate the initial
254
 
          // value of the F_integral above and we therefore just grab
255
 
          // it from the restricted coeffcients
256
 
          u(row) = _ufcs[stage][0]->w()[_coefficient_index[stage][0]][local_to_local_dofs[row]];
257
 
        }
258
 
        t_impl_update_F.stop();
259
 
 
260
 
        // Start iterations
261
 
        while (!newton_converged && newton_iteration < maxiter)
262
 
        {
263
 
          if (_retabulate_J || !reuse_jacobian)
264
 
          {
265
 
            // Tabulate Jacobian
266
 
            Timer t_impl_tt_J("Implicit stage: tabulate_tensor (J)");
267
 
            J_integral.tabulate_tensor(_ufcs[stage][1]->A.data(),
268
 
                                       _ufcs[stage][1]->w(),
269
 
                                       vertex_coordinates.data(),
270
 
                                       local_vert);
271
 
            t_impl_tt_J.stop();
272
 
 
273
 
            // Extract vertex dofs from tabulated tensor
274
 
            Timer t_impl_update_J("Implicit stage: update_J");
275
 
            for (unsigned int row = 0; row < N; row++)
276
 
            {
277
 
              for (unsigned int col=0; col < N; col++)
278
 
              {
279
 
                _J(row, col)
280
 
                  = _ufcs[stage][1]->A[local_to_local_dofs[row]*dof_offset*N
281
 
                                       + local_to_local_dofs[col]];
282
 
              }
283
 
            }
284
 
            t_impl_update_J.stop();
285
 
 
286
 
            // LU factorize Jacobian
287
 
            Timer lu_factorize("Implicit stage: LU factorize");
288
 
            _J_LU.compute(_J);
289
 
            _retabulate_J = false;
290
 
          }
291
 
 
292
 
          // Perform linear solve By forward backward substitution
293
 
          Timer forward_backward_substitution("Implicit stage: fb substitution");
294
 
          dx = _J_LU.solve(F);
295
 
 
296
 
          forward_backward_substitution.stop();
297
 
 
298
 
          // Compute resdiual
299
 
          if (convergence_criterion == "residual")
300
 
            residual = F.norm();
301
 
          else if (convergence_criterion == "incremental")
302
 
            residual = dx.norm();
303
 
          else
304
 
            error("Unknown Newton convergence criterion");
305
 
 
306
 
          // If initial residual
307
 
          if (newton_iteration == 0)
308
 
            residual0 = residual;
309
 
 
310
 
          // Relative residual
311
 
          relative_residual = residual / residual0;
312
 
 
313
 
          // Update solution
314
 
          if (std::abs(1.0 - relaxation) < DOLFIN_EPS)
315
 
            u -= dx;
316
 
          else
317
 
            u -= relaxation*dx;
318
 
 
319
 
          // Update number of iterations
320
 
          ++newton_iteration;
321
 
 
322
 
          // Put solution back into restricted coefficients before
323
 
          // tabulate new residual
324
 
          for (unsigned int row = 0; row < N; row++)
325
 
          {
326
 
            _ufcs[stage][0]->w()[_coefficient_index[stage][0]][local_to_local_dofs[row]]
327
 
              = u(row);
328
 
          }
329
 
          // Tabulate new residual
330
 
          t_impl_tt_F.start();
331
 
          F_integral.tabulate_tensor(_ufcs[stage][0]->A.data(),
332
 
                                     _ufcs[stage][0]->w(),
333
 
                                     vertex_coordinates.data(),
334
 
                                     local_vert);
335
 
          t_impl_tt_F.stop();
336
 
 
337
 
          t_impl_update_F.start();
338
 
          // Extract vertex dofs from tabulated tensor
339
 
          for (unsigned int row = 0; row < N; row++)
340
 
            F(row) = _ufcs[stage][0]->A[local_to_local_dofs[row]];
341
 
          t_impl_update_F.stop();
342
 
 
343
 
          // Output iteration number and residual (only first vertex)
344
 
          if (report && (newton_iteration > 0) && (vert_ind == 0))
345
 
          {
346
 
            info("Point solver newton iteration %d: r (abs) = %.3e (tol = %.3e) "\
347
 
                 "r (rel) = %.3e (tol = %.3e)",
348
 
                 newton_iteration, residual, atol, relative_residual, rtol);
349
 
          }
350
 
 
351
 
          // Check for retabulation of Jacobian
352
 
          if (reuse_jacobian
353
 
              && newton_iteration > iterations_to_retabulate_jacobian
354
 
              &&  !jacobian_retabulated)
355
 
          {
356
 
            jacobian_retabulated = true;
357
 
            _retabulate_J = true;
358
 
 
359
 
            if (vert_ind == 0)
360
 
              info("Retabulating Jacobian.");
361
 
 
362
 
            // If there is a solution coefficient in the jacobian form
363
 
            if (_coefficient_index[stage].size()==2)
364
 
            {
365
 
              // Put solution back into restricted coefficients before
366
 
              // tabulate new jacobian
367
 
              for (unsigned int row=0; row < N; row++)
368
 
                _ufcs[stage][1]->w()[_coefficient_index[stage][1]][local_to_local_dofs[row]] = u(row);
369
 
            }
370
 
 
371
 
          }
372
 
 
373
 
          // Return true if convergence criterion is met
374
 
          if (relative_residual < rtol || residual < atol)
375
 
            newton_converged = true;
376
 
        }
377
 
 
378
 
        if (newton_converged)
379
 
        {
380
 
          Timer t_impl_set("Implicit stage: set");
381
 
          // Put solution back into global stage solution vector
382
 
          _scheme->stage_solutions()[stage]->vector()->set(u.data(), u.size(),
383
 
                                                           &local_to_global_dofs[0]);
384
 
        }
385
 
        else
386
 
        {
387
 
          info("Last iteration before error %d: r (abs) = %.3e (tol = %.3e) "
388
 
               "r (rel) = %.3e (tol = %.3e)", newton_iteration, residual, atol,
389
 
               relative_residual, rtol);
390
 
          error("Newton solver in PointIntegralSolver did not converge.");
391
 
        }
 
177
        _solve_implicit_stage(vert_ind, stage, cell, ufc_cell, vertex_coordinates);
392
178
      }
 
179
 
393
180
    }
394
181
 
395
 
    Timer t_vert_axpy("Step: AXPY solution");
396
 
 
397
 
    // Get local u0 solution and add the stage derivatives
398
 
    _scheme->solution()->vector()->get_local(u0.data(), u0.size(),
399
 
                                             &local_to_global_dofs[0]);
400
 
 
401
 
    // Do the last stage and put back into solution vector
402
 
    FunctionAXPY last_stage = _scheme->last_stage()*dt;
403
 
 
404
 
    // Axpy local solution vectors
405
 
    for (unsigned int stage=0; stage < num_stages; stage++)
406
 
      u0 += last_stage.pairs()[stage].first*local_stage_solutions[stage];
 
182
    Timer t_last_stage("Last stage: tabulate_tensor");
 
183
 
 
184
    // Last stage point integral
 
185
    const ufc::point_integral& integral = *_last_stage_ufc->default_point_integral;
 
186
 
 
187
    // Update coeffcients for last stage
 
188
    // TODO: Pass suitable bool vector here to avoid tabulating all coefficient dofs:
 
189
    _last_stage_ufc->update(cell, vertex_coordinates, ufc_cell);
 
190
                            //integral.enabled_coefficients());
 
191
 
 
192
    // Tabulate cell tensor
 
193
    integral.tabulate_tensor(_last_stage_ufc->A.data(), _last_stage_ufc->w(),
 
194
                             vertex_coordinates.data(),
 
195
                             _vertex_map[vert_ind].second,
 
196
                             ufc_cell.orientation);
 
197
 
 
198
    // Update solution with a tabulation of the last stage
 
199
    for (unsigned int row=0; row < _system_size; row++)
 
200
      _y[row] = _last_stage_ufc->A[_local_to_local_dofs[row]];
407
201
 
408
202
    // Update global solution with last stage
409
 
    _scheme->solution()->vector()->set(u0.data(), local_to_global_dofs.size(),
410
 
                                       &local_to_global_dofs[0]);
411
 
 
412
 
    //p++;
 
203
    _scheme->solution()->vector()->set(_y.data(), _local_to_global_dofs.size(),
 
204
                                       _local_to_global_dofs.data());
413
205
  }
414
206
 
 
207
  for (unsigned int stage=0; stage<_num_stages; stage++)
 
208
    _scheme->stage_solutions()[stage]->vector()->apply("insert");
 
209
 
 
210
  _scheme->solution()->vector()->apply("insert");
 
211
 
415
212
  // Update time
416
213
  *_scheme->t() = t0 + dt;
 
214
 
 
215
}
 
216
//-----------------------------------------------------------------------------
 
217
void PointIntegralSolver::_solve_explicit_stage(std::size_t vert_ind,
 
218
                                                unsigned int stage,
 
219
                                                const ufc::cell& ufc_cell,
 
220
                                                const std::vector<double>& vertex_coordinates)
 
221
{
 
222
  Timer t_expl("Explicit stage");
 
223
 
 
224
  // Local vertex ind
 
225
  const unsigned int local_vert = _vertex_map[vert_ind].second;
 
226
 
 
227
  // Point integral
 
228
  const ufc::point_integral& integral = *_ufcs[stage][0]->default_point_integral;
 
229
 
 
230
  // Tabulate cell tensor
 
231
  Timer t_expl_tt("Explicit stage: tabulate_tensor");
 
232
  integral.tabulate_tensor(_ufcs[stage][0]->A.data(), _ufcs[stage][0]->w(),
 
233
                           vertex_coordinates.data(), local_vert,
 
234
                           ufc_cell.orientation);
 
235
  t_expl_tt.stop();
 
236
 
 
237
  // Extract vertex dofs from tabulated tensor and put them into the local
 
238
  // stage solution vector
 
239
  // Extract vertex dofs from tabulated tensor
 
240
  for (unsigned int row=0; row < _system_size; row++)
 
241
  {
 
242
    _local_stage_solutions[stage][row] = _ufcs[stage][0]->A[_local_to_local_dofs[row]];
 
243
  }
 
244
 
 
245
  // Put solution back into global stage solution vector
 
246
  // NOTE: This so an UFC.update (coefficient restriction) would just work
 
247
  _scheme->stage_solutions()[stage]->vector()->set(
 
248
                 _local_stage_solutions[stage].data(), _system_size,
 
249
                 _local_to_global_dofs.data());
 
250
 
 
251
}
 
252
//-----------------------------------------------------------------------------
 
253
void PointIntegralSolver::_solve_implicit_stage(std::size_t vert_ind,
 
254
                                                unsigned int stage, const Cell& cell,
 
255
                                                const ufc::cell& ufc_cell,
 
256
                                                const std::vector<double>& vertex_coordinates)
 
257
{
 
258
 
 
259
  Timer t_impl("Implicit stage");
 
260
 
 
261
  // Do a simplified newton solve
 
262
  _simplified_newton_solve(vert_ind, stage, cell, ufc_cell, vertex_coordinates);
 
263
 
 
264
  // Put solution back into global stage solution vector
 
265
  _scheme->stage_solutions()[stage]->vector()->set(_local_stage_solutions[stage].data(), \
 
266
                                                   _system_size,
 
267
                                                   _local_to_global_dofs.data());
 
268
 
417
269
}
418
270
//-----------------------------------------------------------------------------
419
271
void PointIntegralSolver::step_interval(double t0, double t1, double dt)
448
300
  }
449
301
}
450
302
//-----------------------------------------------------------------------------
 
303
void PointIntegralSolver::_compute_jacobian(std::vector<double>& jac,
 
304
                                            const std::vector<double>& u,
 
305
                                            unsigned int local_vert,
 
306
                                            UFC& loc_ufc, const Cell& cell,
 
307
                                            const ufc::cell& ufc_cell,
 
308
                                            int coefficient_index,
 
309
                                            const std::vector<double>& vertex_coordinates)
 
310
{
 
311
  const ufc::point_integral& J_integral = *loc_ufc.default_point_integral;
 
312
 
 
313
  //Timer _timer_compute_jac("Implicit stage: Compute jacobian");
 
314
  //Timer t_impl_update("Update_cell");
 
315
  // TODO: Pass suitable bool vector here to avoid tabulating all coefficient dofs:
 
316
  loc_ufc.update(cell, vertex_coordinates, ufc_cell);
 
317
                 //J_integral.enabled_coefficients());
 
318
  //t_impl_update.stop();
 
319
 
 
320
  // If there is a solution coefficient in the jacobian form
 
321
  if (coefficient_index>0)
 
322
  {
 
323
    // Put solution back into restricted coefficients before tabulate new jacobian
 
324
    for (unsigned int row=0; row < _system_size; row++)
 
325
    {
 
326
      loc_ufc.w()[coefficient_index][_local_to_local_dofs[row]] = u[row];
 
327
    }
 
328
  }
 
329
 
 
330
  // Tabulate Jacobian
 
331
  Timer t_impl_tt_jac("Implicit stage: tabulate_tensor (J)");
 
332
  J_integral.tabulate_tensor(loc_ufc.A.data(), loc_ufc.w(),
 
333
                             vertex_coordinates.data(),
 
334
                             local_vert,
 
335
                             ufc_cell.orientation);
 
336
  t_impl_tt_jac.stop();
 
337
 
 
338
  // Extract vertex dofs from tabulated tensor
 
339
  //Timer t_impl_update_jac("Implicit stage: update_jac");
 
340
  for (unsigned int row=0; row < _system_size; row++)
 
341
    for (unsigned int col=0; col < _system_size; col++)
 
342
      jac[row*_system_size + col] = loc_ufc.A[_local_to_local_dofs[row]*
 
343
                                              _dof_offset*_system_size +
 
344
                                              _local_to_local_dofs[col]];
 
345
  //t_impl_update_jac.stop();
 
346
 
 
347
  // LU factorize Jacobian
 
348
  //Timer lu_factorize("Implicit stage: LU factorize");
 
349
  _lu_factorize(jac);
 
350
  _num_jacobian_computations += 1;
 
351
}
 
352
//-----------------------------------------------------------------------------
 
353
void PointIntegralSolver::_lu_factorize(std::vector<double>& A)
 
354
{
 
355
 
 
356
  // Local variables
 
357
  double sum;
 
358
  const int system_size = _system_size;
 
359
 
 
360
  for (int k = 1; k < system_size; k++)
 
361
  {
 
362
 
 
363
    for (int i = 0; i <= k-1; ++i)
 
364
    {
 
365
 
 
366
      sum = 0.0;
 
367
      for (int r = 0; r <= i-1; r++)
 
368
      {
 
369
        sum += A[i*_system_size+r]*A[r*_system_size+k];
 
370
      }
 
371
 
 
372
      A[i*_system_size+k] -=sum;
 
373
 
 
374
      sum = 0.0;
 
375
      for (int r = 0; r <= i-1; r++)
 
376
      {
 
377
        sum += A[k*_system_size+r]*A[r*_system_size+i];
 
378
      }
 
379
 
 
380
      A[k*_system_size+i] = (A[k*_system_size+i]-sum)/A[i*_system_size+i];
 
381
 
 
382
    }
 
383
 
 
384
    sum = 0.0;
 
385
    for (int r = 0; r <= k-1; r++)
 
386
    {
 
387
      sum += A[k*_system_size+r]*A[r*_system_size+k];
 
388
    }
 
389
 
 
390
    A[k*_system_size+k] -= sum;
 
391
 
 
392
  }
 
393
}
 
394
//-----------------------------------------------------------------------------
 
395
void PointIntegralSolver::_forward_backward_subst(const std::vector<double>& A,
 
396
                                                  const std::vector<double>& b,
 
397
                                                  std::vector<double>& x) const
 
398
{
 
399
  // solves Ax = b with forward backward substitution, provided that
 
400
  // A is already LU factorized
 
401
 
 
402
  double sum;
 
403
 
 
404
  x[0] = b[0];
 
405
 
 
406
  // Forward
 
407
  for (unsigned int i=1; i < _system_size; ++i)
 
408
  {
 
409
    sum = 0.0;
 
410
    for (unsigned int j=0; j <= i-1; ++j)
 
411
    {
 
412
      sum = sum + A[i*_system_size+j]*x[j];
 
413
    }
 
414
 
 
415
    x[i] = b[i] -sum;
 
416
  }
 
417
 
 
418
  const unsigned int _system_size_m_1 = _system_size-1;
 
419
  x[_system_size_m_1] = x[_system_size_m_1]/\
 
420
    A[_system_size_m_1*_system_size+_system_size_m_1];
 
421
 
 
422
  // Backward
 
423
  for (int i = _system_size-2; i >= 0; i--)
 
424
  {
 
425
    sum = 0;
 
426
    for (unsigned int j=i+1; j < _system_size; ++j)
 
427
    {
 
428
      sum = sum +A[i*_system_size+j]*x[j];
 
429
    }
 
430
 
 
431
    x[i] = (x[i]-sum)/A[i*_system_size+i];
 
432
  }
 
433
}
 
434
//-----------------------------------------------------------------------------
 
435
double PointIntegralSolver::_norm(const std::vector<double>& vec) const
 
436
{
 
437
  double l2_norm = 0;
 
438
 
 
439
  for (unsigned int i = 0; i < vec.size(); ++i)
 
440
    l2_norm += vec[i]*vec[i];
 
441
 
 
442
  l2_norm = std::sqrt(l2_norm);
 
443
  return l2_norm;
 
444
}
 
445
//-----------------------------------------------------------------------------
451
446
void PointIntegralSolver::_check_forms()
452
447
{
453
448
  // Iterate over stage forms and check they include point integrals
454
 
  std::vector<std::vector<boost::shared_ptr<const Form> > >& stage_forms
 
449
  std::vector<std::vector<std::shared_ptr<const Form> > >& stage_forms
455
450
    = _scheme->stage_forms();
456
451
  for (unsigned int i=0; i < stage_forms.size(); i++)
457
452
  {
466
461
      }
467
462
 
468
463
      // Num dofs per vertex
469
 
      const Mesh& mesh = *stage_forms[i][j]->function_space(0)->mesh();
470
 
      const GenericDofMap& dofmap
471
 
        = *stage_forms[i][j]->function_space(0)->dofmap();
 
464
      const GenericDofMap& dofmap = *stage_forms[i][j]->function_space(0)->dofmap();
472
465
      const unsigned int dofs_per_vertex = dofmap.num_entity_dofs(0);
473
 
      const unsigned int vert_per_cell
474
 
        = mesh.topology()(mesh.topology().dim(), 0).size(0);
 
466
      const unsigned int vert_per_cell = _mesh.topology()(_mesh.topology().dim(), 0).size(0);
475
467
 
476
468
      if (vert_per_cell*dofs_per_vertex != dofmap.max_cell_dimension())
477
469
      {
486
478
void PointIntegralSolver::_init()
487
479
{
488
480
  // Get stage forms
489
 
  std::vector<std::vector<boost::shared_ptr<const Form> > >& stage_forms
 
481
  std::vector<std::vector<std::shared_ptr<const Form> > >& stage_forms
490
482
    = _scheme->stage_forms();
491
483
 
 
484
  // Init local stage solutions
 
485
  for (unsigned int stage=0; stage < _num_stages; stage++)
 
486
    _local_stage_solutions[stage].resize(_system_size);
 
487
 
492
488
  // Init coefficient index and ufcs
493
489
  _coefficient_index.resize(stage_forms.size());
494
490
  _ufcs.resize(stage_forms.size());
496
492
  // Initiate jacobian matrices
497
493
  if (_scheme->implicit())
498
494
  {
499
 
    const unsigned int N
500
 
      = stage_forms[0][0]->function_space(0)->dofmap()->num_entity_dofs(0);
501
 
    _J.resize(N, N);
 
495
    // First count the number of distinct jacobians
 
496
    int max_jacobian_index = 0;
 
497
    for (unsigned int stage=0; stage < _num_stages; stage++)
 
498
      max_jacobian_index = std::max(_scheme->jacobian_index(stage),
 
499
                                    max_jacobian_index);
 
500
 
 
501
    // Create memory for jacobians
 
502
    _jacobians.resize(max_jacobian_index+1);
 
503
    for (int i=0; i<=max_jacobian_index; i++)
 
504
      _jacobians[i].resize(_system_size*_system_size);
 
505
    _recompute_jacobian.resize(max_jacobian_index+1, true);
 
506
 
502
507
  }
503
508
 
 
509
  // Create last stage UFC form
 
510
  _last_stage_ufc = std::make_shared<UFC>(*_scheme->last_stage());
 
511
 
504
512
  // Iterate over stages and collect information
505
513
  for (unsigned int stage = 0; stage < stage_forms.size(); stage++)
506
514
  {
507
515
    // Create a UFC object for first form
508
 
    _ufcs[stage].push_back(boost::make_shared<UFC>(*stage_forms[stage][0]));
 
516
    _ufcs[stage].push_back(std::make_shared<UFC>(*stage_forms[stage][0]));
509
517
 
510
518
    //  If implicit stage
511
519
    if (stage_forms[stage].size()==2)
512
520
    {
513
521
      // Create a UFC object for second form
514
 
      _ufcs[stage].push_back(boost::make_shared<UFC>(*stage_forms[stage][1]));
 
522
      _ufcs[stage].push_back(std::make_shared<UFC>(*stage_forms[stage][1]));
515
523
 
516
524
      // Find coefficient index for each of the two implicit forms
517
525
      for (unsigned int i = 0; i < 2; i++)
529
537
          }
530
538
        }
531
539
      }
 
540
 
 
541
      // Check that nonlinear form includes a coefficient index
 
542
      // (otherwise it might be some error in the form)
 
543
      if (_coefficient_index[stage].size() == 0)
 
544
      {
 
545
        dolfin_error("PointIntegralSolver.cpp",
 
546
                     "constructing PointIntegralSolver",
 
547
                     "Expecting nonlinear form of stage: %d to be dependent "\
 
548
                     "on the stage solution.", stage);
 
549
      }
 
550
 
532
551
    }
533
552
  }
534
553
 
535
 
  // Extract mesh
536
 
  const Mesh& mesh = stage_forms[0][0]->mesh();
537
 
  _vertex_map.resize(mesh.num_vertices());
 
554
  // Build vertex map
 
555
  _vertex_map.resize(_mesh.num_vertices());
538
556
 
539
557
  // Init mesh connections
540
 
  mesh.init(0);
541
 
  const unsigned int dim_t = mesh.topology().dim();
 
558
  _mesh.init(0);
 
559
  const unsigned int dim_t = _mesh.topology().dim();
542
560
 
543
561
  // Iterate over vertices and collect cell and local vertex information
544
 
  for (VertexIterator vert(mesh); !vert.end(); ++vert)
 
562
  for (VertexIterator vert(_mesh); !vert.end(); ++vert)
545
563
  {
546
564
    // First look for cell where the vert is local vert 0
547
565
    bool local_vert_found = false;
561
579
    // local cell 0 and find what local vert the global vert corresponds to
562
580
    if (!local_vert_found)
563
581
    {
564
 
      const Cell cell0(mesh, vert->entities(dim_t)[0]);
 
582
      const Cell cell0(_mesh, vert->entities(dim_t)[0]);
565
583
      _vertex_map[vert->index()].first = cell0.index();
566
584
 
567
585
      unsigned int local_vert_index = 0;
582
600
  }
583
601
}
584
602
//-----------------------------------------------------------------------------
 
603
void PointIntegralSolver::_simplified_newton_solve(
 
604
                              std::size_t vert_ind, unsigned int stage,
 
605
//                            UFC& loc_ufc,
 
606
//                            int coefficient_index_F, int coefficient_index_J,
 
607
                              const Cell& cell,
 
608
                              const ufc::cell& ufc_cell,
 
609
                              const std::vector<double>& vertex_coordinates)
 
610
{
 
611
 
 
612
  //Timer _timer_newton_solve("Implicit stage: Newton solve");
 
613
  const Parameters& newton_solver_params = parameters("newton_solver");
 
614
  const size_t report_vertex = newton_solver_params["report_vertex"];
 
615
  const double kappa = newton_solver_params["kappa"];
 
616
  const double rtol = newton_solver_params["relative_tolerance"];
 
617
  std::size_t max_iterations = newton_solver_params["maximum_iterations"];
 
618
  const double max_relative_previous_residual = newton_solver_params["max_relative_previous_residual"];
 
619
  const double relaxation = newton_solver_params["relaxation_parameter"];
 
620
  const bool report = newton_solver_params["report"];
 
621
  const bool verbose_report = newton_solver_params["verbose_report"];
 
622
  bool always_recompute_jacobian = newton_solver_params["always_recompute_jacobian"];
 
623
  const unsigned int local_vert = _vertex_map[vert_ind].second;
 
624
  UFC& loc_ufc_F = *_ufcs[stage][0];
 
625
  UFC& loc_ufc_J = *_ufcs[stage][1];
 
626
  const int coefficient_index_F = _coefficient_index[stage][0];
 
627
  const int coefficient_index_J = _coefficient_index[stage].size()==2 ?
 
628
    _coefficient_index[stage][1] : -1;
 
629
  const unsigned int jac_index = _scheme->jacobian_index(stage);
 
630
  std::vector<double>& jac = _jacobians[jac_index];
 
631
 
 
632
  if (newton_solver_params["recompute_jacobian_each_solve"])
 
633
    _recompute_jacobian[jac_index] = true;
 
634
 
 
635
  bool newton_solve_restared = false;
 
636
  unsigned int newton_iterations = 0;
 
637
  double relative_previous_residual = 1.0, residual, prev_residual = 1.0,
 
638
    initial_residual = 1.0, relative_residual = 1.0;
 
639
 
 
640
  // Get point integrals
 
641
  const ufc::point_integral& F_integral = *loc_ufc_F.default_point_integral;
 
642
 
 
643
  // Local solution
 
644
  std::vector<double>& u = _local_stage_solutions[stage];
 
645
 
 
646
  // Update with previous local solution and make a backup of solution to
 
647
  // be used in a potential restarting of newton solver
 
648
  for (unsigned int row=0; row < _system_size; row++)
 
649
    _u0[row] = u[row] = loc_ufc_F.w()[coefficient_index_F][_local_to_local_dofs[row]];
 
650
 
 
651
  do
 
652
  {
 
653
 
 
654
    // Tabulate residual
 
655
    Timer t_impl_tt_F("Implicit stage: tabulate_tensor (F)");
 
656
    F_integral.tabulate_tensor(loc_ufc_F.A.data(), loc_ufc_F.w(),
 
657
                               vertex_coordinates.data(),
 
658
                               local_vert,
 
659
                               ufc_cell.orientation);
 
660
    t_impl_tt_F.stop();
 
661
 
 
662
    // Extract vertex dofs from tabulated tensor, together with the old stage
 
663
    // solution
 
664
    for (unsigned int row=0; row < _system_size; row++)
 
665
      _residual[row] = loc_ufc_F.A[_local_to_local_dofs[row]];
 
666
 
 
667
    residual = _norm(_residual);
 
668
    if (newton_iterations == 0)
 
669
      initial_residual = residual;//std::max(residual, DOLFIN_EPS);
 
670
 
 
671
    relative_residual = residual/initial_residual;
 
672
 
 
673
    // Check for relative residual convergence
 
674
    if (relative_residual < rtol)
 
675
      break;
 
676
 
 
677
    // Should we recompute jacobian
 
678
    if (_recompute_jacobian[jac_index] || always_recompute_jacobian)
 
679
    {
 
680
      _compute_jacobian(jac, u, local_vert, loc_ufc_J, cell, ufc_cell,
 
681
                        coefficient_index_J, vertex_coordinates);
 
682
      _recompute_jacobian[jac_index] = false;
 
683
    }
 
684
 
 
685
    // Perform linear solve By forward backward substitution
 
686
    //Timer forward_backward_substitution("Implicit stage: fb substituion");
 
687
    _forward_backward_subst(jac, _residual, _dx);
 
688
    //forward_backward_substitution.stop();
 
689
 
 
690
    // Newton_Iterations == 0
 
691
    if (newton_iterations == 0)
 
692
    {
 
693
      // On first iteration we need an approximation of eta. We take
 
694
      // the one from previous step and increase it slightly. This is
 
695
      // important for linear problems which only should recquire 1
 
696
      // iteration to converge.
 
697
      _eta = _eta > DOLFIN_EPS ? _eta : DOLFIN_EPS;
 
698
      _eta = std::pow(_eta, 0.8);
 
699
    }
 
700
 
 
701
    // 2nd time around
 
702
    else
 
703
    {
 
704
 
 
705
      // How fast are we converging?
 
706
      relative_previous_residual = residual/prev_residual;
 
707
 
 
708
      if (always_recompute_jacobian)
 
709
      {
 
710
 
 
711
        if ((report && vert_ind == report_vertex) || verbose_report)
 
712
          info("Newton solver after %d iterations. vertex: %d, "        \
 
713
               "relative_previous_residual: %.3f, "                     \
 
714
               "relative_residual: %.3e, residual: %.3e.",
 
715
               newton_iterations, vert_ind, relative_previous_residual,
 
716
               relative_residual, residual);
 
717
 
 
718
      }
 
719
 
 
720
      // If we diverge
 
721
      else if (relative_previous_residual > 1)
 
722
      {
 
723
 
 
724
        if ((report && vert_ind == report_vertex) || verbose_report)
 
725
          info("Newton solver diverges after %d iterations. vertex: %d, "               \
 
726
               "relative_previous_residual: %.3f, "                     \
 
727
               "relative_residual: %.3e, residual: %.3e.",
 
728
               newton_iterations, vert_ind, relative_previous_residual,
 
729
               relative_residual, residual);
 
730
 
 
731
        // If we have not restarted newton solve previously
 
732
        if (!newton_solve_restared)
 
733
        {
 
734
          if ((report && vert_ind == report_vertex) || verbose_report)
 
735
            info("Restarting newton solve for vertex: %d", vert_ind);
 
736
 
 
737
          // Reset flags
 
738
          newton_solve_restared = true;
 
739
          always_recompute_jacobian = true;
 
740
 
 
741
          // Reset solution
 
742
          for (unsigned int row=0; row < _system_size; row++)
 
743
            loc_ufc_F.w()[coefficient_index_F][_local_to_local_dofs[row]] = u[row] = _u0[row];
 
744
 
 
745
          // Update variables
 
746
          _eta = parameters("newton_solver")["eta_0"];
 
747
          newton_iterations = 0;
 
748
          relative_previous_residual = prev_residual = initial_residual = relative_residual = 1.0;
 
749
          max_iterations = 400;
 
750
          continue;
 
751
        }
 
752
 
 
753
      }
 
754
 
 
755
      // We converge too slow
 
756
      else if (relative_previous_residual >= max_relative_previous_residual ||
 
757
               residual > (kappa*rtol*(1 - relative_previous_residual)/ \
 
758
                           std::pow(relative_previous_residual,         \
 
759
                                    max_iterations - newton_iterations)))
 
760
      {
 
761
 
 
762
        if ((report && vert_ind == report_vertex) || verbose_report)
 
763
          info("Newton solver converges too slow at iteration %d. vertex: %d, " \
 
764
               "relative_previous_residual: %.3f, "                     \
 
765
               "relative_residual: %.3e, residual: %.3e. Recomputing jacobian.",
 
766
               newton_iterations, vert_ind, relative_previous_residual,
 
767
               relative_residual, residual);
 
768
 
 
769
        _recompute_jacobian[jac_index] = true;
 
770
 
 
771
      }
 
772
      else
 
773
      {
 
774
        if ((report && vert_ind == report_vertex) || verbose_report)
 
775
          info("Newton solver after %d iterations. vertex: %d, "        \
 
776
               "relative_previous_residual: %.3f, "                     \
 
777
               "relative_residual: %.3e, residual: %.3e.",
 
778
               newton_iterations, vert_ind, relative_previous_residual,
 
779
               relative_residual, residual);
 
780
 
 
781
        // Update eta
 
782
        _eta = relative_previous_residual/(1.0 - relative_previous_residual);
 
783
 
 
784
      }
 
785
 
 
786
    }
 
787
 
 
788
    // No convergence
 
789
    if (newton_iterations > max_iterations)
 
790
    {
 
791
      if (report)
 
792
        info("Newton solver did not converge after %d iterations. vertex: %d, " \
 
793
             "relative_previous_residual: %.3f, "                       \
 
794
             "relative_residual: %.3e, residual: %.3e.", max_iterations, vert_ind,
 
795
             relative_previous_residual, relative_residual, residual);
 
796
 
 
797
      error("Newton solver in PointIntegralSolver exeeded maximal iterations.");
 
798
    }
 
799
 
 
800
    // Update solution
 
801
    if (std::abs(1.0 - relaxation) < DOLFIN_EPS)
 
802
      for (unsigned int i=0; i < u.size(); i++)
 
803
        u[i] -= _dx[i];
 
804
    else
 
805
      for (unsigned int i=0; i < u.size(); i++)
 
806
        u[i] -= relaxation*_dx[i];
 
807
 
 
808
    // Put solution back into restricted coefficients before tabulate new residual
 
809
    for (unsigned int row=0; row < _system_size; row++)
 
810
      loc_ufc_F.w()[coefficient_index_F][_local_to_local_dofs[row]] = u[row];
 
811
 
 
812
    prev_residual = residual;
 
813
    newton_iterations++;
 
814
 
 
815
  } while(_eta*relative_residual >= kappa*rtol);
 
816
 
 
817
  if ((report && vert_ind == report_vertex) || verbose_report)
 
818
    info("Newton solver converged after %d iterations. vertex: %d, "\
 
819
         "relative_previous_residual: %.3f, relative_residual: %.3e, "\
 
820
         "residual: %.3e.", newton_iterations, vert_ind,
 
821
         relative_previous_residual, relative_residual, residual);
 
822
 
 
823
}
 
824
//-----------------------------------------------------------------------------