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
|
<?xml version="1.0" encoding="UTF-8"?>
<grammar xmlns:a="http://relaxng.org/ns/compatibility/annotations/1.0" xmlns="http://relaxng.org/ns/structure/1.0" datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
<include href="spud_base.rng"/>
<define name="adjoint_options">
<element name="adjoint">
<a:documentation>Turn on the solution of the adjoint problem.
To compute the adjoint you must specify code for the functional
derivative. If the functional is also specified, then it will be
output in the .stat file.</a:documentation>
<oneOrMore>
<element name="functional">
<a:documentation>A functional to be computed.
If the functional_value is supplied by the user, then the functional gets
evaluated and printed in the .stat files.
If the functional_derivative is supplied by the user, then that is always used.
for the right-hand side of the adjoint computation.
If no functional_derivative is supplied, then automatic differentiation is applied
to compute it from the functional_value.
In short, you must supply at least one of functional_value or functional_derivative,
and you may supply both.</a:documentation>
<attribute name="name">
<data type="string"/>
</attribute>
<optional>
<element name="functional_value">
<a:documentation>The value of the functional to be computed.</a:documentation>
<element name="algorithm">
<a:documentation>Python code for the functional. This code is given the following variables:
n, the current timestep number
time, the start of the timestep
dt, the current size of the timestep
states, which represents the states through time, and may be accessed like
states[0]["MaterialPhase"].vector_fields["Coordinate"], or
states[n-1]["MyOtherMaterialPhase"].scalar_fields["Pressure"], etc.
states will only contain those fields that have been specified in the functional_dependencies element.
This code must compute a real number J, which is the value of the functional associated with this particular
timestep.
Here is a simple example for the functional that evaluates the L2-norm of a field "LayerThickness" at time T=1:
<span font_desc="monospace 10" foreground="blue">
J = 0.0
T = 1.0 # the time at which to evaluate
if time &lt; T &lt;= time+dt:
import numpy
eta_prev = states[n-1]["Fluid"].scalar_fields["LayerThickness"]
eta = states[n]["Fluid"].scalar_fields["LayerThickness"]
# We want to temporally interpolate to evaluate eta at t=1.0
alpha = (time + dt - T) / dt
assert 0 &lt;= alpha &lt; 1
tmp_eta = alpha * eta_prev.val + (1-alpha) * eta.val
# Now we want to integrate that over space
coord = states[0]["Fluid"].vector_fields["Coordinate"]
assert eta.element_count == eta_prev.element_count == coord.element_count
for ele in range(coord.element_count):
t = Transform(ele, coord)
shape = eta_prev.ele_shape(ele)
mass = t.shape_shape(shape, shape)
nodes = eta_prev.ele_nodes(ele)
J = J + numpy.dot(tmp_eta[nodes], numpy.dot(mass, tmp_eta[nodes]))
</span>
<span weight="bold">
If you intend to use automatic differentiation, and want to use
primitives such as sin, cos, exp, etc., you must use those
found in the uncertainties.unumpy package, not those
found in the numpy or math packages.
</span></a:documentation>
<attribute name="name">
<value>functional_value</value>
</attribute>
<ref name="python_code"/>
</element>
<element name="reduction">
<a:documentation>The functional does not exist at any particular point in time; however, the code above is called at every timestep.
This is to allow for temporal localisation. Suppose you wished to compute an integral in time: this allows you to perform
the integral over each element of your time domain separately, and then they must be combined to compute the full value of
the functional. In this example, the relevant reduction would be a sum.</a:documentation>
<element name="sum">
<a:documentation>Sum all of the functional components associated with each timestep to compute the functional</a:documentation>
<empty/>
</element>
</element>
</element>
</optional>
<optional>
<element name="functional_derivative">
<a:documentation>The derivative of the functional at each time level.</a:documentation>
<element name="algorithm">
<a:documentation>Python code for the functional derivative. This code is given the following variables:
n, the timelevel of the field to differentiate with respect to
times, an array containing the value of time for each timelevel
states, which represents the states through time, and may be accessed like
states[0]["MaterialPhase"].vector_fields["Coordinate"], or
states[n]["MyOtherMaterialPhase"].scalar_fields["Pressure"], etc.
states will only contain those fields that have been specified in the functional_dependencies element.
derivative, which is a scalar/vector/tensor field associated with the derivative to be computed.
This code must set the entries of derivative.val. The code should check which variable we are differentiating
with respect to by inspecting derivative.name.</a:documentation>
<attribute name="name">
<value>functional_derivative</value>
</attribute>
<ref name="python_code"/>
</element>
</element>
</optional>
<element name="functional_dependencies">
<a:documentation>The dependencies of the functional at each time level.</a:documentation>
<element name="algorithm">
<a:documentation> Python code for the functional dependencies.
This code defines a function that informs the model what variables at which time levels
will be necessary for the functional computation at this point.
Here is a simple example for the functional that evaluates the L2-norm of a field "LayerThickness" at time T=1:
<span font_desc="monospace 10" foreground="blue">
def dependencies(times, timestep):
if times[0] &lt; 1.0 &lt;= times[1]:
return {"Fluid::Coordinate": [0],
"Fluid::LayerThickness": [timestep-1, timestep]}
else:
return {}
</span>
</a:documentation>
<attribute name="name">
<value>functional_dependencies</value>
</attribute>
<ref name="python_code"/>
</element>
</element>
<optional>
<element name="disable_adjoint_run">
<a:documentation>Disables the adjoint run for this functional.
The functional values are still computed during the forward run.
The resulting stat file is equivalent to the one with this option disabled,
but all adjoint related entries will be zeroed.</a:documentation>
<empty/>
</element>
</optional>
</element>
</oneOrMore>
<optional>
<element name="controls">
<a:documentation>Specifies any control variables </a:documentation>
<zeroOrMore>
<element name="control">
<attribute name="name">
<data type="string"/>
</attribute>
<choice>
<element name="type">
<a:documentation>Initial condition </a:documentation>
<attribute name="name">
<value>initial_condition</value>
</attribute>
<attribute name="field_name">
<data type="string"/>
</attribute>
<ref name="comment"/>
</element>
<element name="type">
<a:documentation>Source term. Note that the source term has to exist in both the forward and in the adjoint run (using the exist_in_both flag)</a:documentation>
<attribute name="name">
<value>source_term</value>
</attribute>
<attribute name="field_name">
<data type="string"/>
</attribute>
<ref name="comment"/>
</element>
</choice>
<optional>
<element name="bounds">
<a:documentation>Defines bounds on the control variables.</a:documentation>
<optional>
<element name="upper_bound">
<a:documentation>Defines an upper bound for this control. The field name must be a a valid field in the state with the same type (e.g. scalar, vector or tensor) and mesh than the control.</a:documentation>
<attribute name="field_name">
<data type="string"/>
</attribute>
<ref name="comment"/>
</element>
</optional>
<optional>
<element name="lower_bound">
<a:documentation>Defines an lower bound for this control. The field name must be a a valid field in the state with the same type (e.g. scalar, vector or tensor) and mesh than the control.</a:documentation>
<attribute name="field_name">
<data type="string"/>
</attribute>
<ref name="comment"/>
</element>
</optional>
</element>
</optional>
</element>
</zeroOrMore>
<optional>
<element name="load_controls">
<a:documentation>Load controls values from the control files.
This is used for optimality only and should be left disabled.</a:documentation>
<empty/>
</element>
</optional>
</element>
</optional>
<optional>
<element name="debug">
<a:documentation>Debugging options for adjoint model development</a:documentation>
<optional>
<element name="replay_forward_run">
<a:documentation>Rerun the forward equation using libadjoint.
This is a debugging option to check if the libadjoint callbacks are implemented correctly.</a:documentation>
<empty/>
</element>
</optional>
<optional>
<element name="html_output">
<a:documentation>Activate the visualisation of the forward and adjoint equations solved.
This option creates an html file for both the forward and adjoint system.</a:documentation>
<empty/>
</element>
</optional>
<optional>
<element name="check_action_transposes">
<a:documentation>Check the transposes of all action callbacks.</a:documentation>
<empty/>
</element>
</optional>
<optional>
<element name="check_action_derivative">
<a:documentation>Check all supplied derivative callbacks of nonlinear operators.</a:documentation>
<empty/>
</element>
</optional>
</element>
</optional>
</element>
</define>
<define name="adjoint_storage">
<element name="adjoint_storage">
<a:documentation>Informs the model whether the field should be computed only in the
forward model, the adjoint model, or in both.
If /adjoint is not enabled, this has no effect whatsoever.</a:documentation>
<choice>
<element name="exists_in_forward">
<a:documentation>Include the field only in the forward state.</a:documentation>
<empty/>
</element>
<element name="exists_in_both">
<a:documentation>Include the field in both the forward and adjoint states.</a:documentation>
<empty/>
</element>
<element name="exists_in_adjoint">
<a:documentation>Include the field only in the adjoint state.</a:documentation>
<empty/>
</element>
</choice>
</element>
</define>
</grammar>
|