~fluidity-core/fluidity/adjoint

« back to all changes in this revision

Viewing changes to tests/shallow_water_optimisation_check_gradient_revolve_2d/adjoint_A.swml

Merge the fluidity_revolve branch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
<?xml version='1.0' encoding='utf-8'?>
 
2
<shallow_water_options>
 
3
  <simulation_name>
 
4
    <string_value lines="1">wave_A</string_value>
 
5
  </simulation_name>
 
6
  <geometry>
 
7
    <dimension>
 
8
      <integer_value rank="0">3</integer_value>
 
9
    </dimension>
 
10
    <mesh name="CoordinateMesh">
 
11
      <from_file file_name="src/mesh_A">
 
12
        <format name="triangle"/>
 
13
        <dimension>
 
14
          <integer_value rank="0">2</integer_value>
 
15
        </dimension>
 
16
        <stat>
 
17
          <include_in_stat/>
 
18
        </stat>
 
19
      </from_file>
 
20
    </mesh>
 
21
    <mesh name="VelocityMesh">
 
22
      <from_mesh>
 
23
        <mesh name="PeriodicMesh"/>
 
24
        <mesh_continuity>
 
25
          <string_value>discontinuous</string_value>
 
26
        </mesh_continuity>
 
27
        <stat>
 
28
          <exclude_from_stat/>
 
29
        </stat>
 
30
      </from_mesh>
 
31
    </mesh>
 
32
    <mesh name="PressureMesh">
 
33
      <from_mesh>
 
34
        <mesh name="PeriodicMesh"/>
 
35
        <mesh_shape>
 
36
          <polynomial_degree>
 
37
            <integer_value rank="0">2</integer_value>
 
38
          </polynomial_degree>
 
39
        </mesh_shape>
 
40
        <stat>
 
41
          <exclude_from_stat/>
 
42
        </stat>
 
43
      </from_mesh>
 
44
    </mesh>
 
45
    <mesh name="PeriodicMesh">
 
46
      <from_mesh>
 
47
        <mesh name="CoordinateMesh"/>
 
48
        <periodic_boundary_conditions name="periodicity">
 
49
          <physical_boundary_ids>
 
50
            <integer_value shape="1" rank="1">10</integer_value>
 
51
          </physical_boundary_ids>
 
52
          <aliased_boundary_ids>
 
53
            <integer_value shape="1" rank="1">8</integer_value>
 
54
          </aliased_boundary_ids>
 
55
          <coordinate_map>
 
56
            <string_value lines="20" type="python">def val(X,t):
 
57
  result = list(X)
 
58
  result[0]=result[0]-1.0
 
59
  return result</string_value>
 
60
          </coordinate_map>
 
61
        </periodic_boundary_conditions>
 
62
        <periodic_boundary_conditions name="periodic_y">
 
63
          <physical_boundary_ids>
 
64
            <integer_value shape="1" rank="1">7</integer_value>
 
65
          </physical_boundary_ids>
 
66
          <aliased_boundary_ids>
 
67
            <integer_value shape="1" rank="1">9</integer_value>
 
68
          </aliased_boundary_ids>
 
69
          <coordinate_map>
 
70
            <string_value lines="20" type="python">def val(X,t):
 
71
  result = list(X)
 
72
  result[1]=result[1]-1.0
 
73
  return result</string_value>
 
74
          </coordinate_map>
 
75
        </periodic_boundary_conditions>
 
76
        <stat>
 
77
          <exclude_from_stat/>
 
78
        </stat>
 
79
      </from_mesh>
 
80
    </mesh>
 
81
    <quadrature>
 
82
      <degree>
 
83
        <integer_value rank="0">4</integer_value>
 
84
      </degree>
 
85
    </quadrature>
 
86
  </geometry>
 
87
  <io>
 
88
    <dump_format>
 
89
      <string_value>vtk</string_value>
 
90
    </dump_format>
 
91
    <dump_period_in_timesteps>
 
92
      <constant>
 
93
        <integer_value rank="0">1</integer_value>
 
94
      </constant>
 
95
    </dump_period_in_timesteps>
 
96
    <output_mesh name="CoordinateMesh"/>
 
97
  </io>
 
98
  <timestepping>
 
99
    <current_time>
 
100
      <real_value rank="0">0</real_value>
 
101
    </current_time>
 
102
    <timestep>
 
103
      <real_value rank="0">0.25</real_value>
 
104
    </timestep>
 
105
    <nonlinear_iterations>
 
106
      <integer_value rank="0">1</integer_value>
 
107
    </nonlinear_iterations>
 
108
    <finish_time>
 
109
      <real_value rank="0">1</real_value>
 
110
    </finish_time>
 
111
  </timestepping>
 
112
  <physical_parameters>
 
113
    <gravity>
 
114
      <magnitude>
 
115
        <real_value rank="0">1</real_value>
 
116
      </magnitude>
 
117
      <vector_field name="GravityDirection" rank="1">
 
118
        <prescribed>
 
119
          <mesh name="CoordinateMesh"/>
 
120
          <value name="WholeMesh">
 
121
            <constant>
 
122
              <real_value shape="3" dim1="dim" rank="1">0 0 -1</real_value>
 
123
            </constant>
 
124
          </value>
 
125
          <output/>
 
126
          <stat>
 
127
            <include_in_stat/>
 
128
          </stat>
 
129
          <detectors>
 
130
            <exclude_from_detectors/>
 
131
          </detectors>
 
132
          <adjoint_storage>
 
133
            <exists_in_forward/>
 
134
          </adjoint_storage>
 
135
        </prescribed>
 
136
      </vector_field>
 
137
    </gravity>
 
138
  </physical_parameters>
 
139
  <material_phase name="Fluid">
 
140
    <vector_field name="Velocity" rank="1">
 
141
      <prognostic>
 
142
        <mesh name="VelocityMesh"/>
 
143
        <equation name="ShallowWater"/>
 
144
        <spatial_discretisation>
 
145
          <discontinuous_galerkin>
 
146
            <advection_scheme>
 
147
              <none/>
 
148
            </advection_scheme>
 
149
          </discontinuous_galerkin>
 
150
          <conservative_advection>
 
151
            <real_value rank="0">0</real_value>
 
152
          </conservative_advection>
 
153
        </spatial_discretisation>
 
154
        <solver>
 
155
          <iterative_method name="gmres">
 
156
            <restart>
 
157
              <integer_value rank="0">500</integer_value>
 
158
            </restart>
 
159
          </iterative_method>
 
160
          <preconditioner name="sor"/>
 
161
          <relative_error>
 
162
            <real_value rank="0">1e-12</real_value>
 
163
          </relative_error>
 
164
          <max_iterations>
 
165
            <integer_value rank="0">500</integer_value>
 
166
          </max_iterations>
 
167
          <never_ignore_solver_failures/>
 
168
          <diagnostics>
 
169
            <monitors/>
 
170
          </diagnostics>
 
171
        </solver>
 
172
        <initial_condition name="WholeMesh">
 
173
          <python>
 
174
            <string_value lines="20" type="python">def val(X, t):
 
175
  import constants_revolve_2d
 
176
  return constants_revolve_2d.u_exact(X, t)</string_value>
 
177
          </python>
 
178
        </initial_condition>
 
179
        <output/>
 
180
        <stat>
 
181
          <include_in_stat/>
 
182
          <previous_time_step>
 
183
            <exclude_from_stat/>
 
184
          </previous_time_step>
 
185
          <nonlinear_field>
 
186
            <exclude_from_stat/>
 
187
          </nonlinear_field>
 
188
        </stat>
 
189
        <convergence>
 
190
          <include_in_convergence/>
 
191
        </convergence>
 
192
        <detectors>
 
193
          <include_in_detectors/>
 
194
        </detectors>
 
195
        <steady_state>
 
196
          <include_in_steady_state/>
 
197
        </steady_state>
 
198
        <consistent_interpolation/>
 
199
      </prognostic>
 
200
    </vector_field>
 
201
    <scalar_field name="LayerThickness" rank="0">
 
202
      <prognostic>
 
203
        <mesh name="PressureMesh"/>
 
204
        <spatial_discretisation>
 
205
          <continuous_galerkin>
 
206
            <advection_terms>
 
207
              <exclude_advection_terms/>
 
208
            </advection_terms>
 
209
          </continuous_galerkin>
 
210
          <conservative_advection>
 
211
            <real_value rank="0">0</real_value>
 
212
          </conservative_advection>
 
213
        </spatial_discretisation>
 
214
        <temporal_discretisation>
 
215
          <theta>
 
216
            <real_value rank="0">0.5</real_value>
 
217
          </theta>
 
218
          <relaxation>
 
219
            <real_value rank="0">1</real_value>
 
220
          </relaxation>
 
221
        </temporal_discretisation>
 
222
        <solver>
 
223
          <iterative_method name="cg"/>
 
224
          <preconditioner name="sor"/>
 
225
          <relative_error>
 
226
            <real_value rank="0">1e-12</real_value>
 
227
          </relative_error>
 
228
          <max_iterations>
 
229
            <integer_value rank="0">500</integer_value>
 
230
          </max_iterations>
 
231
          <never_ignore_solver_failures/>
 
232
          <cache_solver_context/>
 
233
          <diagnostics>
 
234
            <monitors/>
 
235
          </diagnostics>
 
236
        </solver>
 
237
        <initial_condition name="WholeMesh">
 
238
          <python>
 
239
            <string_value lines="20" type="python">def val(X, t):
 
240
  import constants_revolve_2d
 
241
  #return constants_revolve_2d.eta_exact(X, t)
 
242
  return 0.0</string_value>
 
243
          </python>
 
244
        </initial_condition>
 
245
        <mean_layer_thickness>
 
246
          <real_value rank="0">1</real_value>
 
247
        </mean_layer_thickness>
 
248
        <output/>
 
249
        <stat/>
 
250
        <consistent_interpolation/>
 
251
      </prognostic>
 
252
    </scalar_field>
 
253
    <scalar_field name="AnalyticalLayerThickness" rank="0">
 
254
      <prescribed>
 
255
        <mesh name="PressureMesh"/>
 
256
        <value name="WholeMesh">
 
257
          <python>
 
258
            <string_value lines="20" type="python">def val(X, t):
 
259
  import constants_revolve_2d
 
260
  return constants_revolve_2d.eta_exact(X, t)</string_value>
 
261
          </python>
 
262
        </value>
 
263
        <output/>
 
264
        <stat/>
 
265
        <detectors>
 
266
          <exclude_from_detectors/>
 
267
        </detectors>
 
268
        <adjoint_storage>
 
269
          <exists_in_forward/>
 
270
        </adjoint_storage>
 
271
      </prescribed>
 
272
    </scalar_field>
 
273
    <scalar_field name="LayerThicknessError" rank="0">
 
274
      <diagnostic>
 
275
        <algorithm source_field_2_type="scalar" name="scalar_difference" source_field_1_name="LayerThickness" source_field_2_name="AnalyticalLayerThickness" material_phase_support="single" source_field_1_type="scalar">
 
276
          <absolute_difference/>
 
277
        </algorithm>
 
278
        <mesh name="PressureMesh"/>
 
279
        <output/>
 
280
        <stat/>
 
281
        <convergence>
 
282
          <include_in_convergence/>
 
283
        </convergence>
 
284
        <detectors>
 
285
          <include_in_detectors/>
 
286
        </detectors>
 
287
        <steady_state>
 
288
          <include_in_steady_state/>
 
289
        </steady_state>
 
290
        <adjoint_storage>
 
291
          <exists_in_forward/>
 
292
        </adjoint_storage>
 
293
      </diagnostic>
 
294
    </scalar_field>
 
295
    <vector_field name="AnalyticalVelocity" rank="1">
 
296
      <prescribed>
 
297
        <mesh name="VelocityMesh"/>
 
298
        <value name="WholeMesh">
 
299
          <python>
 
300
            <string_value lines="20" type="python">def val(X, t):
 
301
  import constants_revolve_2d
 
302
  return constants_revolve_2d.u_exact(X, t)</string_value>
 
303
          </python>
 
304
        </value>
 
305
        <output/>
 
306
        <stat>
 
307
          <include_in_stat/>
 
308
        </stat>
 
309
        <detectors>
 
310
          <exclude_from_detectors/>
 
311
        </detectors>
 
312
        <adjoint_storage>
 
313
          <exists_in_forward/>
 
314
        </adjoint_storage>
 
315
      </prescribed>
 
316
    </vector_field>
 
317
    <vector_field name="VelocityError" rank="1">
 
318
      <diagnostic>
 
319
        <algorithm source_field_2_type="vector" name="vector_difference" source_field_1_name="Velocity" source_field_2_name="AnalyticalVelocity" material_phase_support="single" source_field_1_type="vector">
 
320
          <absolute_difference/>
 
321
        </algorithm>
 
322
        <mesh name="VelocityMesh"/>
 
323
        <output/>
 
324
        <stat>
 
325
          <include_in_stat/>
 
326
        </stat>
 
327
        <convergence>
 
328
          <include_in_convergence/>
 
329
        </convergence>
 
330
        <detectors>
 
331
          <include_in_detectors/>
 
332
        </detectors>
 
333
        <steady_state>
 
334
          <include_in_steady_state/>
 
335
        </steady_state>
 
336
        <adjoint_storage>
 
337
          <exists_in_forward/>
 
338
        </adjoint_storage>
 
339
      </diagnostic>
 
340
    </vector_field>
 
341
  </material_phase>
 
342
  <adjoint>
 
343
    <checkpointing>
 
344
      <max_memory_checkpoints>
 
345
        <integer_value rank="0">3</integer_value>
 
346
      </max_memory_checkpoints>
 
347
      <max_disk_checkpoints>
 
348
        <integer_value rank="0">0</integer_value>
 
349
      </max_disk_checkpoints>
 
350
    </checkpointing>
 
351
    <functional name="integral_eta_t1">
 
352
      <functional_value>
 
353
        <algorithm name="functional_value">
 
354
          <string_value lines="20" type="python">J = 0.0
 
355
T1 = 0.5 # the time at which to evaluate
 
356
T2 = 1.0 # the time at which to evaluate
 
357
T = -1
 
358
if time &lt; T1 &lt;= time+dt:
 
359
  T = T1
 
360
elif time &lt; T2 &lt;= time+dt:
 
361
  T = T2
 
362
if time &lt; T &lt;= time+dt:
 
363
  import numpy
 
364
  eta_prev = states[n-1]["Fluid"].scalar_fields["LayerThickness"]
 
365
  eta      = states[n]["Fluid"].scalar_fields["LayerThickness"]
 
366
  eta_exact_prev = states[n-1]["Fluid"].scalar_fields["AnalyticalLayerThickness"]
 
367
  eta_exact      = states[n]["Fluid"].scalar_fields["AnalyticalLayerThickness"]
 
368
    
 
369
  # We want to temporally interpolate to evaluate eta at t=1.0
 
370
  alpha = (time + dt - T) / dt
 
371
  assert 0 &lt;= alpha &lt; 1
 
372
  tmp_eta = alpha * (eta_prev.val - eta_exact_prev.val) + (1-alpha) * (eta.val - eta_exact.val)
 
373
  
 
374
  # Now we want to integrate that over space
 
375
  coord = states[0]["Fluid"].vector_fields["Coordinate"]
 
376
  assert eta.element_count == eta_prev.element_count == coord.element_count
 
377
  for ele in range(coord.element_count):
 
378
    t = Transform(ele, coord)
 
379
    shape = eta_prev.ele_shape(ele)
 
380
    mass = t.shape_shape(shape, shape)
 
381
    nodes = eta_prev.ele_nodes(ele)
 
382
    J = J + numpy.dot(tmp_eta[nodes], numpy.dot(mass, tmp_eta[nodes]))</string_value>
 
383
        </algorithm>
 
384
        <reduction>
 
385
          <sum/>
 
386
        </reduction>
 
387
      </functional_value>
 
388
      <functional_dependencies>
 
389
        <algorithm name="functional_dependencies">
 
390
          <string_value lines="20" type="python">def dependencies(times, timestep):
 
391
  if times[0] &lt; 0.5 &lt;= times[1] or times[0] &lt; 1.0 &lt;= times[1]:
 
392
    return {"Fluid::Coordinate": [0],
 
393
            "Fluid::LayerThickness": [timestep-1, timestep],
 
394
            "Fluid::AnalyticalLayerThickness": [timestep-1, timestep]
 
395
            }
 
396
  else:
 
397
    return {}</string_value>
 
398
        </algorithm>
 
399
      </functional_dependencies>
 
400
    </functional>
 
401
    <controls>
 
402
      <control name="InitVel">
 
403
        <type field_name="Fluid::Velocity" name="initial_condition"/>
 
404
      </control>
 
405
      <load_controls/>
 
406
    </controls>
 
407
  </adjoint>
 
408
</shallow_water_options>