1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
// Copyright (C) 2005-2006 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
// First added: 2005-01-27
// Last changed: 2006-07-05
#ifndef __MULTI_ADAPTIVE_TIME_SLAB_H
#define __MULTI_ADAPTIVE_TIME_SLAB_H
#include <dolfin/log/dolfin_log.h>
#include <dolfin/common/types.h>
#include <dolfin/la/uBlasVector.h>
#include "Alloc.h"
#include "Partition.h"
#include "MultiAdaptivity.h"
#include "TimeSlab.h"
namespace dolfin
{
class ODE;
class Method;
class TimeSlabSolver;
/// This class represents a multi-adaptive time slab holding the
/// degrees of freedom for the solution of an ODE between two
/// synchronized time levels a and b, with individual time steps for
/// the different components of the system.
class MultiAdaptiveTimeSlab : public TimeSlab
{
public:
/// Constructor
MultiAdaptiveTimeSlab(ODE& ode);
/// Destructor
~MultiAdaptiveTimeSlab();
/// Build time slab, return end time
real build(real a, real b);
/// Solve time slab system
bool solve();
/// Check if current solution can be accepted
bool check(bool first);
/// Shift time slab (prepare for next time slab)
bool shift();
/// Reset to initial data
void reset();
/// Prepare sample at time t
void sample(real t);
/// Sample solution value of given component at given time
real usample(uint i, real t);
/// Sample time step size for given component at given time
real ksample(uint i, real t);
/// Sample residual for given component at given time
real rsample(uint i, real t);
/// Display time slab data
void disp() const;
/// Friends
friend class MultiAdaptiveFixedPointSolver;
friend class MultiAdaptiveNewtonSolver;
friend class MultiAdaptiveJacobian;
friend class UpdatedMultiAdaptiveJacobian;
friend class MultiAdaptivePreconditioner;
friend class MultiAdaptivity;
private:
// Reallocate all data
void allocData(real a, real b);
// Create time slab
real createTimeSlab(real a, real b, uint offset);
// Create time slab data
void create_s(real t0, real t1, uint offset, uint end);
void create_e(uint index, uint subslab, real a, real b);
void create_j(uint index);
void create_d(uint index, uint element, uint subslab, real a0, real b0);
// Reallocation of data
void alloc_s(uint newsize);
void alloc_e(uint newsize);
void alloc_j(uint newsize);
void alloc_d(uint newsize);
// Compute length of time slab
real computeEndTime(real a, real b, uint offset, uint& end);
// Compute size of data
real computeDataSize(real a, real b, uint offset);
// Compute number of dependencies to components with smaller time steps
uint countDependencies(uint i0);
// Compute number of dependencies to components with smaller time steps
uint countDependencies(uint i0, real b0);
// Check if the given time is within the given interval
bool within(real t, real a, real b) const;
// Check if the first given interval is within the second interval
bool within(real a0, real b0, real a1, real b1) const;
// Cover all elements in sub slab and return e1, with e0 <= e < e1
uint coverSlab(int subslab, uint e0);
// Cover all elements in next sub slab and return next sub slab
uint coverNext(int subslab, uint element);
// Cover given time for all components
void coverTime(real t);
// Compute maximum of all element residuals
real computeMaxResiduals();
// Evaluate right-hand side at quadrature points of given element (cG)
void cGfeval(real* f, uint s0, uint e0, uint i0, real a0, real b0, real k0);
// Evaluate right-hand side at quadrature points of given element (dG)
void dGfeval(real* f, uint s0, uint e0, uint i0, real a0, real b0, real k0);
// Choose solver
TimeSlabSolver* chooseSolver();
//--- Time slab data ---
real* sa; // Mapping s --> start time t of sub slab s
real* sb; // Mapping s --> end time t of sub slab s
uint* ei; // Mapping e --> component index i of element e
uint* es; // Mapping e --> time slab s containing element e
uint* ee; // Mapping e --> previous element e of element e
uint* ed; // Mapping e --> first dependency d of element e
real* jx; // Mapping j --> value of dof j
int* de; // Mapping d --> element e of dependency d
//--- Size of time slab data ---
Alloc size_s; // Allocation data for sub slabs s
Alloc size_e; // Allocation data for elements e
Alloc size_j; // Allocation data for dofs j
Alloc size_d; // Allocation data for dependencies d
uint ns; // Number of sub slabs
uint ne; // Number of elements
uint nj; // Number of dofs
uint nd; // Number of dependencies
//--- Auxiliary data, size N ---
TimeSlabSolver* solver; // The solver (size N if diagonally damped)
MultiAdaptivity adaptivity; // Adaptive time step regulation (size 3N)
Partition partition; // Time step partitioning (size N)
int* elast; // Last element for each component (size N)
real* f0; // Right-hand side at left end-point for cG (size N)
//--- Auxiliary data ---
uint emax; // Last covered element for sample
real kmin; // Minimum time step (exluding threshold modified)
uBlasVector u; // The interpolated solution at a given time
};
}
#endif
|