forked from pymor/pymor
-
Notifications
You must be signed in to change notification settings - Fork 0
/
basic.py
283 lines (243 loc) · 12.3 KB
/
basic.py
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
import numpy as np
from pymor.algorithms.timestepping import TimeStepper
from pymor.models.interface import Model
from pymor.operators.block import BlockOperatorBase
from pymor.operators.constructions import IdentityOperator, VectorOperator, ZeroOperator
from pymor.vectorarrays.interface import VectorArray
from pymor.vectorarrays.numpy import NumpyVectorSpace
class StationaryModel(Model):
"""Generic class for models of stationary problems.
This class describes discrete problems given by the equation::
L(u(μ), μ) = F(μ)
with a vector-like right-hand side F and a (possibly non-linear) operator L.
Note that even when solving a variational formulation where F is a
functional and not a vector, F has to be specified as a vector-like
|Operator| (mapping scalars to vectors). This ensures that in the complex
case both L and F are anti-linear in the test variable.
Parameters
----------
operator
The |Operator| L.
rhs
The vector F. Either a |VectorArray| of length 1 or a vector-like
|Operator|.
output_functional
|Operator| mapping a given solution to the model output. In many applications,
this will be a |Functional|, i.e. an |Operator| mapping to scalars.
This is not required, however.
products
A dict of inner product |Operators| defined on the discrete space the
problem is posed on. For each product with key `'x'` a corresponding
attribute `x_product`, as well as a norm method `x_norm` is added to
the model.
error_estimator
An error estimator for the problem. This can be any object with
an `estimate_error(U, mu, m)` method. If `error_estimator` is
not `None`, an `estimate_error(U, mu)` method is added to the
model which will call `error_estimator.estimate_error(U, mu, self)`.
visualizer
A visualizer for the problem. This can be any object with
a `visualize(U, m, ...)` method. If `visualizer`
is not `None`, a `visualize(U, *args, **kwargs)` method is added
to the model which forwards its arguments to the
visualizer's `visualize` method.
name
Name of the model.
"""
def __init__(self, operator, rhs, output_functional=None, products=None,
error_estimator=None, visualizer=None, name=None):
if isinstance(rhs, VectorArray):
assert rhs in operator.range
rhs = VectorOperator(rhs, name='rhs')
assert rhs.range == operator.range and rhs.source.is_scalar and rhs.linear
assert output_functional is None or output_functional.source == operator.source
super().__init__(products=products, error_estimator=error_estimator, visualizer=visualizer, name=name)
self.__auto_init(locals())
self.solution_space = operator.source
self.linear = operator.linear and (output_functional is None or output_functional.linear)
if output_functional is not None:
self.dim_output = output_functional.range.dim
def __str__(self):
return (
f'{self.name}\n'
f' class: {self.__class__.__name__}\n'
f' {"linear" if self.linear else "non-linear"}\n'
f' solution_space: {self.solution_space}\n'
f' dim_output: {self.dim_output}\n'
)
def _compute_solution(self, mu=None, **kwargs):
return self.operator.apply_inverse(self.rhs.as_range_array(mu), mu=mu)
def _compute_solution_d_mu_single_direction(self, parameter, index, solution, mu):
lhs_d_mu = self.operator.d_mu(parameter, index).apply(solution, mu=mu)
rhs_d_mu = self.rhs.d_mu(parameter, index).as_range_array(mu)
rhs = rhs_d_mu - lhs_d_mu
return self.operator.jacobian(solution, mu=mu).apply_inverse(rhs)
_compute_allowed_kwargs = frozenset({'use_adjoint'})
def _compute_output_d_mu(self, solution, mu, return_array=False, use_adjoint=None):
"""compute the gradient of the output functional w.r.t. the parameters
Parameters
----------
solution
Internal model state for the given |Parameter value|
mu
|Parameter value| for which to compute the gradient
return_array
if `True`, return the output gradient as a |NumPy array|.
Otherwise, return a dict of gradients for each |Parameter|.
use_adjoint
if `None` use standard approach, if `True`, use
the adjoint solution for a more efficient way of computing the gradient.
See Section 1.6.2 in [HPUU09]_ for more details.
So far, the adjoint approach is only valid for linear models.
Returns
-------
The gradient as a |NumPy array| or a dict of |NumPy arrays|.
"""
if use_adjoint is None:
use_adjoint = True if (self.output_functional.linear and self.operator.linear) else False
if not use_adjoint:
return super()._compute_output_d_mu(solution, mu, return_array)
else:
assert self.output_functional is not None
assert self.operator.linear
assert self.output_functional.linear
dual_solutions = self.operator.range.empty()
for d in range(self.output_functional.range.dim):
dual_problem = self.with_(operator=self.operator.H, rhs=self.output_functional.H.as_range_array(mu)[d])
dual_solutions.append(dual_problem.solve(mu))
gradients = [] if return_array else {}
for (parameter, size) in self.parameters.items():
array = np.empty(shape=(size,self.output_functional.range.dim))
for index in range(size):
output_partial_dmu = self.output_functional.d_mu(parameter, index).apply(solution,
mu=mu).to_numpy()[0]
lhs_d_mu = self.operator.d_mu(parameter, index).apply2(dual_solutions, solution, mu=mu)[:,0]
rhs_d_mu = self.rhs.d_mu(parameter, index).apply_adjoint(dual_solutions, mu=mu).to_numpy()[:,0]
array[index] = output_partial_dmu + rhs_d_mu - lhs_d_mu
if return_array:
gradients.extend(array)
else:
gradients[parameter] = array
if return_array:
return np.array(gradients)
else:
return gradients
class InstationaryModel(Model):
"""Generic class for models of instationary problems.
This class describes instationary problems given by the equations::
M * ∂_t u(t, μ) + L(u(μ), t, μ) = F(t, μ)
u(0, μ) = u_0(μ)
for t in [0,T], where L is a (possibly non-linear) time-dependent
|Operator|, F is a time-dependent vector-like |Operator|, and u_0 the
initial data. The mass |Operator| M is assumed to be linear.
Parameters
----------
T
The final time T.
initial_data
The initial data `u_0`. Either a |VectorArray| of length 1 or
(for the |Parameter|-dependent case) a vector-like |Operator|
(i.e. a linear |Operator| with `source.dim == 1`) which
applied to `NumpyVectorArray(np.array([1]))` will yield the
initial data for given |parameter values|.
operator
The |Operator| L.
rhs
The right-hand side F.
mass
The mass |Operator| `M`. If `None`, the identity is assumed.
time_stepper
The :class:`time-stepper <pymor.algorithms.timestepping.TimeStepper>`
to be used by :meth:`~pymor.models.interface.Model.solve`.
num_values
The number of returned vectors of the solution trajectory. If `None`, each
intermediate vector that is calculated is returned.
output_functional
|Operator| mapping a given solution to the model output. In many applications,
this will be a |Functional|, i.e. an |Operator| mapping to scalars.
This is not required, however.
products
A dict of product |Operators| defined on the discrete space the
problem is posed on. For each product with key `'x'` a corresponding
attribute `x_product`, as well as a norm method `x_norm` is added to
the model.
error_estimator
An error estimator for the problem. This can be any object with
an `estimate_error(U, mu, m)` method. If `error_estimator` is
not `None`, an `estimate_error(U, mu)` method is added to the
model which will call `error_estimator.estimate_error(U, mu, self)`.
visualizer
A visualizer for the problem. This can be any object with
a `visualize(U, m, ...)` method. If `visualizer`
is not `None`, a `visualize(U, *args, **kwargs)` method is added
to the model which forwards its arguments to the
visualizer's `visualize` method.
name
Name of the model.
"""
def __init__(self, T, initial_data, operator, rhs, mass=None, time_stepper=None, num_values=None,
output_functional=None, products=None, error_estimator=None, visualizer=None, name=None):
if isinstance(rhs, VectorArray):
assert rhs in operator.range
rhs = VectorOperator(rhs, name='rhs')
if isinstance(initial_data, VectorArray):
assert initial_data in operator.source
initial_data = VectorOperator(initial_data, name='initial_data')
mass = mass or IdentityOperator(operator.source)
rhs = rhs or ZeroOperator(operator.source, NumpyVectorSpace(1))
assert isinstance(time_stepper, TimeStepper)
assert initial_data.source.is_scalar
assert operator.source == initial_data.range
assert rhs.linear and rhs.range == operator.range and rhs.source.is_scalar
assert mass.linear and mass.source == mass.range == operator.source
assert output_functional is None or output_functional.source == operator.source
super().__init__(products=products, error_estimator=error_estimator, visualizer=visualizer, name=name)
self.parameters_internal = {'t': 1}
self.__auto_init(locals())
self.solution_space = operator.source
self.linear = operator.linear and (output_functional is None or output_functional.linear)
if output_functional is not None:
self.dim_output = output_functional.range.dim
def __str__(self):
return (
f'{self.name}\n'
f' class: {self.__class__.__name__}\n'
f' {"linear" if self.linear else "non-linear"}\n'
f' T: {self.T}\n'
f' solution_space: {self.solution_space}\n'
f' dim_output: {self.dim_output}\n'
)
def with_time_stepper(self, **kwargs):
return self.with_(time_stepper=self.time_stepper.with_(**kwargs))
def _compute_solution(self, mu=None, **kwargs):
mu = mu.with_(t=0.)
U0 = self.initial_data.as_range_array(mu)
U = self.time_stepper.solve(operator=self.operator,
rhs=None if isinstance(self.rhs, ZeroOperator) else self.rhs,
initial_data=U0,
mass=None if isinstance(self.mass, IdentityOperator) else self.mass,
initial_time=0, end_time=self.T, mu=mu, num_values=self.num_values)
return U
def to_lti(self):
"""Convert model to |LTIModel|.
This method interprets the given model as an |LTIModel|
in the following way::
- self.operator -> A
self.rhs -> B
self.output_functional -> C
None -> D
self.mass -> E
"""
if self.output_functional is None:
raise ValueError('No output defined.')
A = - self.operator
B = self.rhs
C = self.output_functional
E = self.mass
if not all(op.linear for op in [A, B, C, E]):
raise ValueError('Operators not linear.')
from pymor.models.iosys import LTIModel
return LTIModel(A, B, C, E=E, visualizer=self.visualizer)