@@ -25,20 +25,36 @@ def local_vector(u, *, readonly=False):
25
25
26
26
27
27
class L2Cholesky :
28
+ """Mass matrix Cholesky factorization for a (real) DG space.
29
+
30
+ Parameters
31
+ ----------
32
+
33
+ space : WithGeometry
34
+ DG space.
35
+ constant_jacobian : bool
36
+ Whether the mass matrix is constant.
37
+ """
38
+
28
39
def __init__ (self , space , * , constant_jacobian = True ):
40
+ if fd .utils .complex_mode :
41
+ raise NotImplementedError ("complex not supported" )
42
+
29
43
self ._space = space
30
44
self ._constant_jacobian = constant_jacobian
31
- self ._pc = None
45
+ self ._cached_pc = None
32
46
33
47
@property
34
- def space (self ):
48
+ def space (self ) -> fd .functionspaceimpl .WithGeometry :
49
+ """Function space.
50
+ """
51
+
35
52
return self ._space
36
53
37
- @property
38
- def pc (self ):
54
+ def _pc (self ):
39
55
import petsc4py .PETSc as PETSc
40
56
41
- if self ._pc is None :
57
+ if self ._cached_pc is None :
42
58
M = fd .assemble (fd .inner (fd .TrialFunction (self .space ), fd .TestFunction (self .space )) * fd .dx ,
43
59
mat_type = "aij" )
44
60
M_local = M .petscmat .getDiagonalBlock ()
@@ -50,24 +66,70 @@ def pc(self):
50
66
pc .setUp ()
51
67
52
68
if self ._constant_jacobian :
53
- self ._pc = M , M_local , pc
69
+ self ._cached_pc = M , M_local , pc
54
70
else :
55
- _ , _ , pc = self ._pc
71
+ _ , _ , pc = self ._cached_pc
56
72
57
73
return pc
58
74
59
75
def C_inv_action (self , u ):
76
+ """For the Cholesky factorization
77
+
78
+ ... math :
79
+
80
+ M = C C^T,
81
+
82
+ compute the action of :math:`C^{-1}`.
83
+
84
+ Parameters
85
+ ----------
86
+
87
+ u : Function or Cofunction
88
+ Compute :math:`C^{-1} \t ilde{u}` where :math:`\t ilde{u}` is the
89
+ vector of degrees of freedom for :math:`u`.
90
+
91
+ Returns
92
+ -------
93
+
94
+ v : Cofunction
95
+ Has vector of degrees of freedom :math:`C^{-1} \t ilde{u}`.
96
+ """
97
+
98
+ pc = self ._pc ()
60
99
v = fd .Cofunction (self .space .dual ())
61
100
with u .dat .vec_ro as u_v , v .dat .vec_wo as v_v :
62
101
with local_vector (u_v , readonly = True ) as u_v_s , local_vector (v_v ) as v_v_s :
63
- self . pc .applySymmetricLeft (u_v_s , v_v_s )
102
+ pc .applySymmetricLeft (u_v_s , v_v_s )
64
103
return v
65
104
66
105
def C_T_inv_action (self , u ):
106
+ """For the Cholesky factorization
107
+
108
+ ... math :
109
+
110
+ M = C C^T,
111
+
112
+ compute the action of :math:`C^{-T}`.
113
+
114
+ Parameters
115
+ ----------
116
+
117
+ u : Function or Cofunction
118
+ Compute :math:`C^{-T} \t ilde{u}` where :math:`\t ilde{u}` is the
119
+ vector of degrees of freedom for :math:`u`.
120
+
121
+ Returns
122
+ -------
123
+
124
+ v : Function
125
+ Has vector of degrees of freedom :math:`C^{-T} \t ilde{u}`.
126
+ """
127
+
128
+ pc = self ._pc ()
67
129
v = fd .Function (self .space )
68
130
with u .dat .vec_ro as u_v , v .dat .vec_wo as v_v :
69
131
with local_vector (u_v , readonly = True ) as u_v_s , local_vector (v_v ) as v_v_s :
70
- self . pc .applySymmetricRight (u_v_s , v_v_s )
132
+ pc .applySymmetricRight (u_v_s , v_v_s )
71
133
return v
72
134
73
135
@@ -78,7 +140,7 @@ class L2RieszMap(fd.RieszMap):
78
140
----------
79
141
80
142
target : WithGeometry
81
- Target space.
143
+ Function space.
82
144
83
145
Keyword arguments are passed to the :class:`firedrake.RieszMap`
84
146
constructor.
@@ -91,11 +153,41 @@ def __init__(self, target, **kwargs):
91
153
92
154
93
155
def is_dg_space (space ):
156
+ """Return whether a function space is DG.
157
+
158
+ Parameters
159
+ ----------
160
+
161
+ space : WithGeometry
162
+ The function space.
163
+
164
+ Returns
165
+ -------
166
+
167
+ bool
168
+ Whether the function space is DG.
169
+ """
170
+
94
171
e , _ = finat .element_factory .convert (space .ufl_element ())
95
172
return e .is_dg ()
96
173
97
174
98
175
def dg_space (space ):
176
+ """Construct a DG space containing a given function space as a subspace.
177
+
178
+ Parameters
179
+ ----------
180
+
181
+ space : WithGeometry
182
+ A function space.
183
+
184
+ Returns
185
+ -------
186
+
187
+ WithGeometry
188
+ A DG space containing `space` as a subspace. May be `space`.
189
+ """
190
+
99
191
if is_dg_space (space ):
100
192
return space
101
193
else :
@@ -112,15 +204,13 @@ class L2TransformedFunctional(AbstractReducedFunctional):
112
204
where
113
205
114
206
- :math:`J` is the functional definining an optimization problem.
115
- - :math:`\Pi` is the :math:`L^2` projection from a discontinuous
116
- superspace of the control space.
207
+ - :math:`\Pi` is the :math:`L^2` projection from a DG space containing
208
+ the control space as a subspace .
117
209
- :math:`\Xi` represents a change of basis from an :math:`L^2`
118
- orthonormal basis to the finite element basis for the discontinuous
119
- superspace.
210
+ orthonormal basis to the finite element basis for the DG space.
120
211
121
212
The optimization is therefore transformed into an optimization problem
122
- using an :math:`L^2` orthonormal basis for a discontinuous finite element
123
- space.
213
+ using an :math:`L^2` orthonormal basis for a DG finite element space.
124
214
125
215
The transformation is related to the factorization in section 4.1 of
126
216
https://doi.org/10.1137/18M1175239 -- specifically the factorization
@@ -134,8 +224,8 @@ class L2TransformedFunctional(AbstractReducedFunctional):
134
224
controls : Control or Sequence[Control]
135
225
Controls. Must be :class:`firedrake.Function` objects.
136
226
riesz_map : L2RieszMap or Sequence[L2RieszMap]
137
- Used for projecting from the discontinuous space onto the control
138
- space. Ignored for DG controls.
227
+ Used for projecting from the DG space onto the control space. Ignored
228
+ for DG controls.
139
229
alpha : Real
140
230
Modifies the functional, equivalent to adding an extra term to
141
231
:math:`J \circ \Pi`
@@ -145,8 +235,7 @@ class L2TransformedFunctional(AbstractReducedFunctional):
145
235
\frac{1}{2} \alpha \left\| m_D - \Pi ( m_D ) \right\|_{L^2}^2.
146
236
147
237
e.g. in a minimization problem this adds a penalty term which can
148
- be used to avoid ill-posedness due to the use of a larger discontinuous
149
- space.
238
+ be used to avoid ill-posedness due to the use of a larger DG space.
150
239
tape : Tape
151
240
Tape used in evaluations involving :math:`J`.
152
241
"""
@@ -239,7 +328,7 @@ def map_result(self, m):
239
328
240
329
m : firedrake.Function or Sequence[firedrake.Function]
241
330
The result of the optimization. Represents an expansion in an
242
- :math:`L^2` orthonormal basis for the discontinuous space.
331
+ :math:`L^2` orthonormal basis for the DG space.
243
332
244
333
Returns
245
334
-------
0 commit comments