← Back to branch summary

~fluidity-core/fluidity/sea-ice-branch

« back to all changes in this revision

Viewing changes to tests/shallow_water_adjoint_default_controls_2d/adjoint_D.swml

  • Committer: Simon Mouradian
  • Date: 2012-10-19 10:35:59 UTC
  • mfrom: (3520.32.371 fluidity)
  • Revision ID: simon.mouradian06@imperial.ac.uk-20121019103559-y36qa47phc69q8sc
merge from trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/shallow_water_adjoint_default_controls_2d/adjoint_D.swml
1
 
<?xml version='1.0' encoding='utf-8'?>
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">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="code" language="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="code" language="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>
92
 
      <constant>
93
 
        <real_value rank="0">0.0</real_value>
94
 
      </constant>
95
 
    </dump_period>
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.03125</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">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 0.0 -1.0</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="preonly"/>
156
 
          <preconditioner name="lu"/>
157
 
          <relative_error>
158
 
            <real_value rank="0">1.0e-7</real_value>
159
 
          </relative_error>
160
 
          <max_iterations>
161
 
            <integer_value rank="0">500</integer_value>
162
 
          </max_iterations>
163
 
          <never_ignore_solver_failures/>
164
 
          <diagnostics>
165
 
            <monitors/>
166
 
          </diagnostics>
167
 
        </solver>
168
 
        <initial_condition name="WholeMesh">
169
 
          <python>
170
 
            <string_value lines="20" type="code" language="python">def val(X, t):
171
 
  import constants_2d
172
 
  return constants_2d.u_exact(X, t)</string_value>
173
 
          </python>
174
 
        </initial_condition>
175
 
        <vector_field name="Source" rank="1">
176
 
          <prescribed>
177
 
            <value name="WholeMesh">
178
 
              <python>
179
 
                <string_value lines="20" type="code" language="python">def val(X, t):
180
 
  import constants_2d
181
 
  return constants_2d.u_src(X, t)</string_value>
182
 
              </python>
183
 
            </value>
184
 
            <output/>
185
 
            <stat>
186
 
              <include_in_stat/>
187
 
            </stat>
188
 
            <detectors>
189
 
              <exclude_from_detectors/>
190
 
            </detectors>
191
 
            <adjoint_storage>
192
 
              <exists_in_forward/>
193
 
            </adjoint_storage>
194
 
          </prescribed>
195
 
        </vector_field>
196
 
        <output/>
197
 
        <stat>
198
 
          <include_in_stat/>
199
 
          <previous_time_step>
200
 
            <exclude_from_stat/>
201
 
          </previous_time_step>
202
 
          <nonlinear_field>
203
 
            <exclude_from_stat/>
204
 
          </nonlinear_field>
205
 
        </stat>
206
 
        <convergence>
207
 
          <include_in_convergence/>
208
 
        </convergence>
209
 
        <detectors>
210
 
          <include_in_detectors/>
211
 
        </detectors>
212
 
        <steady_state>
213
 
          <include_in_steady_state/>
214
 
        </steady_state>
215
 
        <consistent_interpolation/>
216
 
      </prognostic>
217
 
    </vector_field>
218
 
    <scalar_field name="LayerThickness" rank="0">
219
 
      <prognostic>
220
 
        <mesh name="PressureMesh"/>
221
 
        <spatial_discretisation>
222
 
          <continuous_galerkin>
223
 
            <advection_terms>
224
 
              <exclude_advection_terms/>
225
 
            </advection_terms>
226
 
          </continuous_galerkin>
227
 
          <conservative_advection>
228
 
            <real_value rank="0">0</real_value>
229
 
          </conservative_advection>
230
 
        </spatial_discretisation>
231
 
        <temporal_discretisation>
232
 
          <theta>
233
 
            <real_value rank="0">0.5</real_value>
234
 
          </theta>
235
 
          <relaxation>
236
 
            <real_value rank="0">1</real_value>
237
 
          </relaxation>
238
 
        </temporal_discretisation>
239
 
        <solver>
240
 
          <iterative_method name="preonly"/>
241
 
          <preconditioner name="lu"/>
242
 
          <relative_error>
243
 
            <real_value rank="0">1.0e-7</real_value>
244
 
          </relative_error>
245
 
          <max_iterations>
246
 
            <integer_value rank="0">500</integer_value>
247
 
          </max_iterations>
248
 
          <never_ignore_solver_failures/>
249
 
          <cache_solver_context/>
250
 
          <diagnostics>
251
 
            <monitors/>
252
 
          </diagnostics>
253
 
        </solver>
254
 
        <initial_condition name="WholeMesh">
255
 
          <python>
256
 
            <string_value lines="20" type="code" language="python">def val(X, t):
257
 
  import constants_2d
258
 
  return constants_2d.eta_exact(X, t)</string_value>
259
 
          </python>
260
 
        </initial_condition>
261
 
        <mean_layer_thickness>
262
 
          <real_value rank="0">0.5</real_value>
263
 
        </mean_layer_thickness>
264
 
        <scalar_field name="Source">
265
 
          <prescribed>
266
 
            <value name="WholeMesh">
267
 
              <python>
268
 
                <string_value lines="20" type="code" language="python">def val(X, t):
269
 
  import constants_2d
270
 
  return constants_2d.eta_src(X, t)</string_value>
271
 
              </python>
272
 
            </value>
273
 
            <output/>
274
 
            <stat/>
275
 
            <detectors>
276
 
              <exclude_from_detectors/>
277
 
            </detectors>
278
 
            <adjoint_storage>
279
 
              <exists_in_forward/>
280
 
            </adjoint_storage>
281
 
          </prescribed>
282
 
        </scalar_field>
283
 
        <output/>
284
 
        <stat/>
285
 
        <consistent_interpolation/>
286
 
      </prognostic>
287
 
    </scalar_field>
288
 
    <scalar_field name="AnalyticalLayerThickness" rank="0">
289
 
      <prescribed>
290
 
        <mesh name="PressureMesh"/>
291
 
        <value name="WholeMesh">
292
 
          <python>
293
 
            <string_value lines="20" type="code" language="python">def val(X, t):
294
 
  import constants_2d
295
 
  return constants_2d.eta_exact(X, t)</string_value>
296
 
          </python>
297
 
        </value>
298
 
        <output/>
299
 
        <stat/>
300
 
        <detectors>
301
 
          <exclude_from_detectors/>
302
 
        </detectors>
303
 
        <adjoint_storage>
304
 
          <exists_in_both/>
305
 
        </adjoint_storage>
306
 
      </prescribed>
307
 
    </scalar_field>
308
 
    <scalar_field name="LayerThicknessError" rank="0">
309
 
      <diagnostic>
310
 
        <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">
311
 
          <absolute_difference/>
312
 
        </algorithm>
313
 
        <mesh name="PressureMesh"/>
314
 
        <output/>
315
 
        <stat/>
316
 
        <convergence>
317
 
          <include_in_convergence/>
318
 
        </convergence>
319
 
        <detectors>
320
 
          <include_in_detectors/>
321
 
        </detectors>
322
 
        <steady_state>
323
 
          <include_in_steady_state/>
324
 
        </steady_state>
325
 
        <adjoint_storage>
326
 
          <exists_in_forward/>
327
 
        </adjoint_storage>
328
 
      </diagnostic>
329
 
    </scalar_field>
330
 
    <scalar_field name="dJdh" rank="0">
331
 
      <diagnostic>
332
 
        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
333
 
          <string_value lines="20" type="code" language="python">import fluidity.state_types
334
 
import numpy
335
 
field.val[:] = 0.0
336
 
dJ = 0.0
337
 
 
338
 
coord  = state.vector_fields["Coordinate"]
339
 
dJdu0 = state.vector_fields['integral_eta_t1_InitVelCtrl_TotalDerivative']
340
 
dJdusrc = state.vector_fields['integral_eta_t1_VelSrcCtrl_TotalDerivative']
341
 
dJdeta0 = state.scalar_fields['integral_eta_t1_InitEtaCtrl_TotalDerivative']
342
 
dJdetasrc = state.scalar_fields['integral_eta_t1_EtaSrcCtrl_TotalDerivative']
343
 
eta_exact = state.scalar_fields['AnalyticalLayerThickness']
344
 
u_exact = state.vector_fields['AnalyticalVelocity']
345
 
eta_src = state.scalar_fields['MyEtaSource']
346
 
u_src = state.vector_fields['MyVelocitySource']
347
 
 
348
 
if time == 0.0:
349
 
  for node in range(dJdeta0.node_count):
350
 
    dJ = dJ + numpy.dot(dJdeta0.node_val(node), eta_exact.node_val(node))
351
 
  for node in range(dJdu0.node_count):    
352
 
    dJ = dJ + numpy.dot(dJdu0.node_val(node), u_exact.node_val(node))
353
 
else:
354
 
  for node in range(dJdetasrc.node_count):
355
 
    dJ = dJ + numpy.dot(dJdetasrc.node_val(node), eta_src.node_val(node))
356
 
  for node in range(dJdusrc.node_count):    
357
 
    dJ = dJ + numpy.dot(dJdusrc.node_val(node), u_src.node_val(node))
358
 
 
359
 
field.val[:] = dJ</string_value>
360
 
          <depends>
361
 
            <string_value lines="1">MyEtaSource,MyVelocitySource</string_value>
362
 
          </depends>
363
 
        </algorithm>
364
 
        <mesh name="VelocityMesh"/>
365
 
        <output/>
366
 
        <stat/>
367
 
        <convergence>
368
 
          <include_in_convergence/>
369
 
        </convergence>
370
 
        <detectors>
371
 
          <include_in_detectors/>
372
 
        </detectors>
373
 
        <steady_state>
374
 
          <include_in_steady_state/>
375
 
        </steady_state>
376
 
        <adjoint_storage>
377
 
          <exists_in_adjoint/>
378
 
        </adjoint_storage>
379
 
      </diagnostic>
380
 
    </scalar_field>
381
 
    <scalar_field name="MyEtaSource" rank="0">
382
 
      <diagnostic>
383
 
        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
384
 
          <string_value lines="20" type="code" language="python">import constants_2d as constants
385
 
timeprime = time + constants.theta*dt
386
 
coord  = state.vector_fields["Coordinate"]
387
 
field[:] = 0.0
388
 
for ele in range(coord.ele_count):
389
 
  ele_vals = [constants.eta_src(x, timeprime) for x in coord.remap_ele(ele, field.mesh)]
390
 
  field.set(field.ele_nodes(ele), ele_vals)</string_value>
391
 
        </algorithm>
392
 
        <mesh name="PressureMesh"/>
393
 
        <output/>
394
 
        <stat/>
395
 
        <convergence>
396
 
          <include_in_convergence/>
397
 
        </convergence>
398
 
        <detectors>
399
 
          <include_in_detectors/>
400
 
        </detectors>
401
 
        <steady_state>
402
 
          <include_in_steady_state/>
403
 
        </steady_state>
404
 
        <adjoint_storage>
405
 
          <exists_in_adjoint/>
406
 
        </adjoint_storage>
407
 
      </diagnostic>
408
 
    </scalar_field>
409
 
    <vector_field name="AnalyticalVelocity" rank="1">
410
 
      <prescribed>
411
 
        <mesh name="VelocityMesh"/>
412
 
        <value name="WholeMesh">
413
 
          <python>
414
 
            <string_value lines="20" type="code" language="python">def val(X, t):
415
 
  import constants_2d
416
 
  return constants_2d.u_exact(X, t)</string_value>
417
 
          </python>
418
 
        </value>
419
 
        <output/>
420
 
        <stat>
421
 
          <include_in_stat/>
422
 
        </stat>
423
 
        <detectors>
424
 
          <exclude_from_detectors/>
425
 
        </detectors>
426
 
        <adjoint_storage>
427
 
          <exists_in_both/>
428
 
        </adjoint_storage>
429
 
      </prescribed>
430
 
    </vector_field>
431
 
    <vector_field name="VelocityError" rank="1">
432
 
      <diagnostic>
433
 
        <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">
434
 
          <absolute_difference/>
435
 
        </algorithm>
436
 
        <mesh name="VelocityMesh"/>
437
 
        <output/>
438
 
        <stat>
439
 
          <include_in_stat/>
440
 
        </stat>
441
 
        <convergence>
442
 
          <include_in_convergence/>
443
 
        </convergence>
444
 
        <detectors>
445
 
          <include_in_detectors/>
446
 
        </detectors>
447
 
        <steady_state>
448
 
          <include_in_steady_state/>
449
 
        </steady_state>
450
 
        <adjoint_storage>
451
 
          <exists_in_forward/>
452
 
        </adjoint_storage>
453
 
      </diagnostic>
454
 
    </vector_field>
455
 
    <vector_field name="MyVelocitySource" rank="1">
456
 
      <diagnostic>
457
 
        <algorithm name="vector_python_diagnostic" material_phase_support="single">
458
 
          <string_value lines="20" type="code" language="python">import constants_2d as constants
459
 
timeprime = time + constants.theta*dt
460
 
coord  = state.vector_fields["Coordinate"]
461
 
field[:] = 0.0
462
 
for ele in range(field.ele_count):
463
 
  ele_vals = [constants.u_src(x, timeprime) for x in coord.remap_ele(ele, field.mesh)]
464
 
  field.set(field.ele_nodes(ele), ele_vals)</string_value>
465
 
        </algorithm>
466
 
        <mesh name="VelocityMesh"/>
467
 
        <output/>
468
 
        <stat>
469
 
          <include_in_stat/>
470
 
        </stat>
471
 
        <convergence>
472
 
          <include_in_convergence/>
473
 
        </convergence>
474
 
        <detectors>
475
 
          <include_in_detectors/>
476
 
        </detectors>
477
 
        <steady_state>
478
 
          <include_in_steady_state/>
479
 
        </steady_state>
480
 
        <adjoint_storage>
481
 
          <exists_in_adjoint/>
482
 
        </adjoint_storage>
483
 
      </diagnostic>
484
 
    </vector_field>
485
 
  </material_phase>
486
 
  <adjoint>
487
 
    <functional name="integral_eta_t1">
488
 
      <functional_value>
489
 
        <algorithm name="functional_value">
490
 
          <string_value lines="20" type="code" language="python">J = 0.0
491
 
T = 1.0 # the time at which to evaluate
492
 
if time &lt; T &lt;= time+dt:
493
 
  import numpy
494
 
  eta_prev = states[n-1]["Fluid"].scalar_fields["LayerThickness"]
495
 
  eta      = states[n]["Fluid"].scalar_fields["LayerThickness"]
496
 
  
497
 
  # We want to temporally interpolate to evaluate eta at t=1.0
498
 
  alpha = (time + dt - T) / dt
499
 
  assert 0 &lt;= alpha &lt; 1
500
 
  tmp_eta = alpha * eta_prev.val + (1-alpha) * eta.val
501
 
  
502
 
  # Now we want to integrate that over space
503
 
  coord = states[0]["Fluid"].vector_fields["Coordinate"]
504
 
  assert eta.element_count == eta_prev.element_count == coord.element_count
505
 
  for ele in range(coord.element_count):
506
 
    t = Transform(ele, coord)
507
 
    shape = eta_prev.ele_shape(ele)
508
 
    mass = t.shape_shape(shape, shape)
509
 
    nodes = eta_prev.ele_nodes(ele)
510
 
    J = J + numpy.dot(tmp_eta[nodes], numpy.dot(mass, tmp_eta[nodes]))</string_value>
511
 
        </algorithm>
512
 
        <reduction>
513
 
          <sum/>
514
 
        </reduction>
515
 
      </functional_value>
516
 
      <functional_dependencies>
517
 
        <algorithm name="functional_dependencies">
518
 
          <string_value lines="20" type="code" language="python">def dependencies(times, timestep):
519
 
  if times[0] &lt; 1.0 &lt;= times[1]:
520
 
    return {"Fluid::Coordinate": [0],
521
 
            "Fluid::LayerThickness": [timestep-1, timestep]}
522
 
  else:
523
 
    return {}</string_value>
524
 
        </algorithm>
525
 
      </functional_dependencies>
526
 
    </functional>
527
 
    <controls>
528
 
      <control name="InitVelCtrl">
529
 
        <type field_name="Velocity" name="initial_condition"/>
530
 
      </control>
531
 
      <control name="InitEtaCtrl">
532
 
        <type field_name="LayerThickness" name="initial_condition"/>
533
 
      </control>
534
 
      <control name="VelSrcCtrl">
535
 
        <type field_name="VelocitySource" name="source_term"/>
536
 
      </control>
537
 
      <control name="EtaSrcCtrl">
538
 
        <type field_name="LayerThicknessSource" name="source_term"/>
539
 
      </control>
540
 
    </controls>
541
 
    <debug>
542
 
      <replay_forward_run/>
543
 
      <html_output/>
544
 
    </debug>
545
 
  </adjoint>
546
 
</shallow_water_options>