~reducedmodelling/fluidity/ROM_Non-intrusive-ann
» Revision
3425.13.1
: tests/shallow_water_optimisation/adjoint_D.swml
« back to all changes in this revision
Viewing changes to tests/shallow_water_optimisation/adjoint_D.swml
- browse files at revision 3425.13.1
- compare with another revision
- download diff
- download tarball
- view history from revision 3425.13.1
- Committer: Simon Funke
- Date: 2011-06-01 19:46:34 UTC
- mto: (3425.14.1 adjoint)
- mto: This revision was merged to the branch mainline in revision 3521.
- Revision ID: simon.funke@gmail.com-20110601194634-pu9tmo4zf6o3btgt
optimisation example and schema's
- files added:
- schemas/optimiser.rnc
- tests/shallow_water_optimisation
- tests/shallow_water_optimisation/src
added
removed
tests/shallow_water_optimisation/adjoint_D.swml
1
<?xml version='1.0' encoding='utf8'?>
2
<shallow_water_options>
3
<simulation_name>
4
<string_value lines="1">wave_D</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_D">
12
<format name="triangle"/>
13
<dimension>
14
<integer_value rank="0">1</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">1</integer_value>
51
</physical_boundary_ids>
52
<aliased_boundary_ids>
53
<integer_value shape="1" rank="1">2</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
<stat>
63
<exclude_from_stat/>
64
</stat>
65
</from_mesh>
66
</mesh>
67
<quadrature>
68
<degree>
69
<integer_value rank="0">4</integer_value>
70
</degree>
71
</quadrature>
72
</geometry>
73
<io>
74
<dump_format>
75
<string_value>vtk</string_value>
76
</dump_format>
77
<dump_period_in_timesteps>
78
<constant>
79
<integer_value rank="0">1</integer_value>
80
</constant>
81
</dump_period_in_timesteps>
82
<output_mesh name="CoordinateMesh"/>
83
</io>
84
<timestepping>
85
<current_time>
86
<real_value rank="0">0</real_value>
87
</current_time>
88
<timestep>
89
<real_value rank="0">0.03125</real_value>
90
</timestep>
91
<nonlinear_iterations>
92
<integer_value rank="0">1</integer_value>
93
</nonlinear_iterations>
94
<finish_time>
95
<real_value rank="0">1.0</real_value>
96
</finish_time>
97
</timestepping>
98
<physical_parameters>
99
<gravity>
100
<magnitude>
101
<real_value rank="0">0.1</real_value>
102
</magnitude>
103
<vector_field name="GravityDirection" rank="1">
104
<prescribed>
105
<mesh name="CoordinateMesh"/>
106
<value name="WholeMesh">
107
<constant>
108
<real_value shape="3" dim1="dim" rank="1">0.0 0.0 -1.0</real_value>
109
</constant>
110
</value>
111
<output/>
112
<stat>
113
<include_in_stat/>
114
</stat>
115
<detectors>
116
<exclude_from_detectors/>
117
</detectors>
118
<adjoint_storage>
119
<exists_in_forward/>
120
</adjoint_storage>
121
</prescribed>
122
</vector_field>
123
</gravity>
124
</physical_parameters>
125
<material_phase name="Fluid">
126
<vector_field name="Velocity" rank="1">
127
<prognostic>
128
<mesh name="VelocityMesh"/>
129
<equation name="ShallowWater"/>
130
<spatial_discretisation>
131
<discontinuous_galerkin>
132
<advection_scheme>
133
<none/>
134
</advection_scheme>
135
</discontinuous_galerkin>
136
<conservative_advection>
137
<real_value rank="0">0</real_value>
138
</conservative_advection>
139
</spatial_discretisation>
140
<solver>
141
<iterative_method name="preonly"/>
142
<preconditioner name="lu"/>
143
<relative_error>
144
<real_value rank="0">1.0e-7</real_value>
145
</relative_error>
146
<max_iterations>
147
<integer_value rank="0">500</integer_value>
148
</max_iterations>
149
<never_ignore_solver_failures/>
150
<diagnostics>
151
<monitors/>
152
</diagnostics>
153
</solver>
154
<initial_condition name="WholeMesh">
155
<python>
156
<string_value lines="20" type="python">def val(X, t):
157
import constants
158
return constants.u_exact(X[0], 0.0)</string_value>
159
</python>
160
</initial_condition>
161
<vector_field name="Source" rank="1">
162
<prescribed>
163
<value name="WholeMesh">
164
<python>
165
<string_value lines="20" type="python">def val(X, t):
166
import constants
167
return constants.u_src(X[0], t)</string_value>
168
</python>
169
</value>
170
<output/>
171
<stat>
172
<include_in_stat/>
173
</stat>
174
<detectors>
175
<exclude_from_detectors/>
176
</detectors>
177
<adjoint_storage>
178
<exists_in_forward/>
179
</adjoint_storage>
180
</prescribed>
181
</vector_field>
182
<output/>
183
<stat>
184
<include_in_stat/>
185
<previous_time_step>
186
<exclude_from_stat/>
187
</previous_time_step>
188
<nonlinear_field>
189
<exclude_from_stat/>
190
</nonlinear_field>
191
</stat>
192
<convergence>
193
<include_in_convergence/>
194
</convergence>
195
<detectors>
196
<include_in_detectors/>
197
</detectors>
198
<steady_state>
199
<include_in_steady_state/>
200
</steady_state>
201
<consistent_interpolation/>
202
</prognostic>
203
</vector_field>
204
<scalar_field name="LayerThickness" rank="0">
205
<prognostic>
206
<mesh name="PressureMesh"/>
207
<spatial_discretisation>
208
<continuous_galerkin>
209
<advection_terms>
210
<exclude_advection_terms/>
211
</advection_terms>
212
</continuous_galerkin>
213
<conservative_advection>
214
<real_value rank="0">0</real_value>
215
</conservative_advection>
216
</spatial_discretisation>
217
<temporal_discretisation>
218
<theta>
219
<real_value rank="0">0.5</real_value>
220
</theta>
221
<relaxation>
222
<real_value rank="0">1</real_value>
223
</relaxation>
224
</temporal_discretisation>
225
<solver>
226
<iterative_method name="preonly"/>
227
<preconditioner name="lu"/>
228
<relative_error>
229
<real_value rank="0">1.0e-7</real_value>
230
</relative_error>
231
<max_iterations>
232
<integer_value rank="0">500</integer_value>
233
</max_iterations>
234
<never_ignore_solver_failures/>
235
<cache_solver_context/>
236
<diagnostics>
237
<monitors/>
238
</diagnostics>
239
</solver>
240
<initial_condition name="WholeMesh">
241
<python>
242
<string_value lines="20" type="python">def val(X, t):
243
import constants
244
return constants.eta_exact(X[0], 0.0)</string_value>
245
</python>
246
</initial_condition>
247
<mean_layer_thickness>
248
<real_value rank="0">0.5</real_value>
249
</mean_layer_thickness>
250
<scalar_field name="Source">
251
<prescribed>
252
<value name="WholeMesh">
253
<python>
254
<string_value lines="20" type="python">def val(X, t):
255
import constants
256
return constants.eta_src(X[0], t)</string_value>
257
</python>
258
</value>
259
<output/>
260
<stat/>
261
<detectors>
262
<exclude_from_detectors/>
263
</detectors>
264
<adjoint_storage>
265
<exists_in_forward/>
266
</adjoint_storage>
267
</prescribed>
268
</scalar_field>
269
<output/>
270
<stat/>
271
<consistent_interpolation/>
272
</prognostic>
273
</scalar_field>
274
<scalar_field name="AnalyticalLayerThickness" rank="0">
275
<prescribed>
276
<mesh name="PressureMesh"/>
277
<value name="WholeMesh">
278
<python>
279
<string_value lines="20" type="python">def val(X,t):
280
import constants
281
x=X[0]
282
return constants.eta_exact(x, t)</string_value>
283
</python>
284
</value>
285
<output/>
286
<stat/>
287
<detectors>
288
<exclude_from_detectors/>
289
</detectors>
290
<adjoint_storage>
291
<exists_in_forward/>
292
</adjoint_storage>
293
</prescribed>
294
</scalar_field>
295
<scalar_field name="LayerThicknessError" rank="0">
296
<diagnostic>
297
<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">
298
<absolute_difference/>
299
</algorithm>
300
<mesh name="PressureMesh"/>
301
<output/>
302
<stat/>
303
<convergence>
304
<include_in_convergence/>
305
</convergence>
306
<detectors>
307
<include_in_detectors/>
308
</detectors>
309
<steady_state>
310
<include_in_steady_state/>
311
</steady_state>
312
<adjoint_storage>
313
<exists_in_forward/>
314
</adjoint_storage>
315
</diagnostic>
316
</scalar_field>
317
<scalar_field name="LayerThicknessSourceDerivative" rank="0">
318
<diagnostic>
319
<algorithm name="scalar_python_diagnostic" material_phase_support="single">
320
<string_value lines="20" type="python">import constants
321
import numpy, numpy.linalg
322
323
coord = state.vector_fields["Coordinate"]
324
eta_mesh = state.meshes["PressureMesh"]
325
u_mesh = state.meshes["VelocityMesh"]
326
timeprime = time + constants.theta*dt
327
field.val[:] = 0.0
328
329
for ele in range(coord.ele_count):
330
t = Transform(ele, coord)
331
shape = field.ele_shape(ele)
332
eta_mass = t.shape_shape(shape, shape)
333
u_shape = u_mesh.shape
334
u_mass = t.shape_shape(u_shape, u_shape)
335
336
if time == 0.0:
337
source = numpy.dot(eta_mass, [constants.eta_exact(x[0], 0.0) for x in coord.remap_ele(ele, field.mesh)])
338
else:
339
# First, the eta source term contribution
340
source = numpy.dot(eta_mass, [constants.eta_src(x[0], timeprime)*abs(dt) for x in coord.remap_ele(ele, field.mesh)])
341
342
# Now, the u source term contribution
343
#Linv_Musrc = numpy.dot(numpy.linalg.inv(u_mass), # + coriolis matrix in here!!!!
344
# numpy.dot(u_mass, [constants.u_src(x, timeprime) for x in coord.remap_ele(ele, u_mesh)]))
345
Linv_Musrc = numpy.array([constants.u_src(x[0], timeprime) for x in coord.remap_ele(ele, u_mesh)])
346
347
u_dshape = t.grad(u_mesh.shape)
348
cT = t.shape_dshape(shape, u_dshape)
349
#print "Initial eta_src: ", source
350
for dim in range(cT.shape[2]):
351
source += numpy.dot(cT[:,:,dim], Linv_Musrc[:,dim]) * (dt**2 * constants.d0 * constants.theta)
352
#print "Final eta_src: ", source
353
field.addto(field.ele_nodes(ele), source)</string_value>
354
</algorithm>
355
<mesh name="PressureMesh"/>
356
<output/>
357
<stat/>
358
<convergence>
359
<include_in_convergence/>
360
</convergence>
361
<detectors>
362
<include_in_detectors/>
363
</detectors>
364
<steady_state>
365
<include_in_steady_state/>
366
</steady_state>
367
<adjoint_storage>
368
<exists_in_adjoint/>
369
</adjoint_storage>
370
</diagnostic>
371
</scalar_field>
372
<scalar_field name="dJdh" rank="0">
373
<diagnostic>
374
<algorithm name="scalar_python_diagnostic" material_phase_support="single">
375
<string_value lines="20" type="python">import fluidity.state_types
376
import numpy
377
378
coord = state.vector_fields["Coordinate"]
379
local_adj_delta_u = state.vector_fields["AdjointLocalVelocityDelta"]
380
du_src = state.vector_fields["VelocitySourceDerivative"]
381
adj_u = state.vector_fields["AdjointVelocity"]
382
383
tmpval = numpy.zeros(du_src.val.shape)
384
adj_delta_u = fluidity.state_types.VectorField("GlobalAdjointVelocityDelta", tmpval, 0, "/none", "None")
385
adj_delta_u.mesh = du_src.mesh
386
for node in range(adj_delta_u.node_count):
387
adj_delta_u.set(node, [-1 * local_adj_delta_u.node_val(node)[0], 0, 0])
388
389
adj_delta_eta = state.scalar_fields["AdjointLayerThicknessDelta"]
390
deta_src = state.scalar_fields["LayerThicknessSourceDerivative"]
391
adj_eta = state.scalar_fields["AdjointLayerThickness"]
392
393
J = 0.0
394
if time == 0.0:
395
J = numpy.dot(adj_eta.val, deta_src.val) + sum([numpy.dot(adj_u.val[:,dim], du_src.val[:,dim]) for dim in range(3)])
396
else:
397
J = numpy.dot(adj_delta_eta.val, deta_src.val) + sum([numpy.dot(adj_delta_u.val[:,dim], du_src.val[:,dim]) for dim in range(3)])
398
399
field.val[:] = J</string_value>
400
<depends>
401
<string_value lines="1">LayerThicknessSourceDerivative,VelocitySourceDerivative</string_value>
402
</depends>
403
</algorithm>
404
<mesh name="VelocityMesh"/>
405
<output/>
406
<stat/>
407
<convergence>
408
<include_in_convergence/>
409
</convergence>
410
<detectors>
411
<include_in_detectors/>
412
</detectors>
413
<steady_state>
414
<include_in_steady_state/>
415
</steady_state>
416
<adjoint_storage>
417
<exists_in_adjoint/>
418
</adjoint_storage>
419
</diagnostic>
420
</scalar_field>
421
<vector_field name="AnalyticalVelocity" rank="1">
422
<prescribed>
423
<mesh name="VelocityMesh"/>
424
<value name="WholeMesh">
425
<python>
426
<string_value lines="20" type="python">def val(X,t):
427
import constants
428
x = X[0]
429
return constants.u_exact(x, t)</string_value>
430
</python>
431
</value>
432
<output/>
433
<stat>
434
<include_in_stat/>
435
</stat>
436
<detectors>
437
<exclude_from_detectors/>
438
</detectors>
439
<adjoint_storage>
440
<exists_in_forward/>
441
</adjoint_storage>
442
</prescribed>
443
</vector_field>
444
<vector_field name="VelocityError" rank="1">
445
<diagnostic>
446
<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">
447
<absolute_difference/>
448
</algorithm>
449
<mesh name="VelocityMesh"/>
450
<output/>
451
<stat>
452
<include_in_stat/>
453
</stat>
454
<convergence>
455
<include_in_convergence/>
456
</convergence>
457
<detectors>
458
<include_in_detectors/>
459
</detectors>
460
<steady_state>
461
<include_in_steady_state/>
462
</steady_state>
463
<adjoint_storage>
464
<exists_in_forward/>
465
</adjoint_storage>
466
</diagnostic>
467
</vector_field>
468
<vector_field name="VelocitySourceDerivative" rank="1">
469
<diagnostic>
470
<algorithm name="vector_python_diagnostic" material_phase_support="single">
471
<string_value lines="20" type="python">import constants
472
import numpy, numpy.linalg
473
474
coord = state.vector_fields["Coordinate"]
475
u_mesh = state.meshes["VelocityMesh"]
476
timeprime = time + constants.theta*dt
477
field.val[:] = 0.0
478
479
for ele in range(coord.ele_count):
480
t = Transform(ele, coord)
481
u_shape = u_mesh.shape
482
u_mass = t.shape_shape(u_shape, u_shape)
483
484
if time == 0.0:
485
source = numpy.dot(u_mass, [constants.u_exact(x[0], 0.0) for x in coord.remap_ele(ele, field.mesh)])
486
else:
487
# Add the u source term contribution
488
source = numpy.dot(u_mass, [constants.u_src(x[0], timeprime)*abs(dt) for x in coord.remap_ele(ele, field.mesh)])
489
#print "Final u_src: ", source
490
field.addto(field.ele_nodes(ele), source)</string_value>
491
</algorithm>
492
<mesh name="VelocityMesh"/>
493
<output/>
494
<stat>
495
<include_in_stat/>
496
</stat>
497
<convergence>
498
<include_in_convergence/>
499
</convergence>
500
<detectors>
501
<include_in_detectors/>
502
</detectors>
503
<steady_state>
504
<include_in_steady_state/>
505
</steady_state>
506
<adjoint_storage>
507
<exists_in_adjoint/>
508
</adjoint_storage>
509
</diagnostic>
510
</vector_field>
511
</material_phase>
512
<adjoint>
513
<functional name="integral_eta_t1">
514
<functional_value>
515
<algorithm name="functional_value">
516
<string_value lines="20" type="python">J = 0.0
517
T = 1.0 # the time at which to evaluate
518
if time < T <= time+dt:
519
import numpy
520
eta_prev = states[n-1]["Fluid"].scalar_fields["LayerThickness"]
521
eta = states[n]["Fluid"].scalar_fields["LayerThickness"]
522
523
# We want to temporally interpolate to evaluate eta at t=1.0
524
alpha = (time + dt - T) / dt
525
assert 0 <= alpha < 1
526
tmp_eta = alpha * eta_prev.val + (1-alpha) * eta.val
527
528
# Now we want to integrate that over space
529
coord = states[0]["Fluid"].vector_fields["Coordinate"]
530
assert eta.element_count == eta_prev.element_count == coord.element_count
531
for ele in range(coord.element_count):
532
t = Transform(ele, coord)
533
shape = eta_prev.ele_shape(ele)
534
mass = t.shape_shape(shape, shape)
535
nodes = eta_prev.ele_nodes(ele)
536
J = J + numpy.dot(tmp_eta[nodes], numpy.dot(mass, tmp_eta[nodes]))</string_value>
537
</algorithm>
538
<reduction>
539
<sum/>
540
</reduction>
541
</functional_value>
542
<functional_dependencies>
543
<algorithm name="functional_dependencies">
544
<string_value lines="20" type="python">def dependencies(times, timestep):
545
if times[0] < 1.0 <= times[1]:
546
return {"Fluid::Coordinate": [0],
547
"Fluid::LayerThickness": [timestep-1, timestep]}
548
else:
549
return {}</string_value>
550
</algorithm>
551
</functional_dependencies>
552
</functional>
553
<debug>
554
<html_output/>
555
</debug>
556
</adjoint>
557
</shallow_water_options>
b'\\ No newline at end of file'
Loggerhead is a web-based interface for Breezy
Version: 2.0.1