← Back to branch summary

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

« back to all changes in this revision

Viewing changes to tests/shallow_water_optimisation/adjoint_C.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_optimisation/adjoint_C.swml
1
 
<?xml version='1.0' encoding='utf-8'?>
2
 
<shallow_water_options>
3
 
  <simulation_name>
4
 
    <string_value lines="1">wave_C</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_C">
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="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
 
        <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.0625</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</real_value>
96
 
    </finish_time>
97
 
  </timestepping>
98
 
  <physical_parameters>
99
 
    <gravity>
100
 
      <magnitude>
101
 
        <real_value rank="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 -1</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">1e-07</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="code" language="python">def val(X, t):
157
 
  import shallow_water_optimisation as constants
158
 
  return constants.u_exact(X[0], 0.0)</string_value>
159
 
          </python>
160
 
        </initial_condition>
161
 
        <output/>
162
 
        <stat>
163
 
          <include_in_stat/>
164
 
          <previous_time_step>
165
 
            <exclude_from_stat/>
166
 
          </previous_time_step>
167
 
          <nonlinear_field>
168
 
            <exclude_from_stat/>
169
 
          </nonlinear_field>
170
 
        </stat>
171
 
        <convergence>
172
 
          <include_in_convergence/>
173
 
        </convergence>
174
 
        <detectors>
175
 
          <include_in_detectors/>
176
 
        </detectors>
177
 
        <steady_state>
178
 
          <include_in_steady_state/>
179
 
        </steady_state>
180
 
        <consistent_interpolation/>
181
 
      </prognostic>
182
 
    </vector_field>
183
 
    <scalar_field name="LayerThickness" rank="0">
184
 
      <prognostic>
185
 
        <mesh name="PressureMesh"/>
186
 
        <spatial_discretisation>
187
 
          <continuous_galerkin>
188
 
            <advection_terms>
189
 
              <exclude_advection_terms/>
190
 
            </advection_terms>
191
 
          </continuous_galerkin>
192
 
          <conservative_advection>
193
 
            <real_value rank="0">0</real_value>
194
 
          </conservative_advection>
195
 
        </spatial_discretisation>
196
 
        <temporal_discretisation>
197
 
          <theta>
198
 
            <real_value rank="0">0.5</real_value>
199
 
          </theta>
200
 
          <relaxation>
201
 
            <real_value rank="0">1</real_value>
202
 
          </relaxation>
203
 
        </temporal_discretisation>
204
 
        <solver>
205
 
          <iterative_method name="preonly"/>
206
 
          <preconditioner name="lu"/>
207
 
          <relative_error>
208
 
            <real_value rank="0">1e-07</real_value>
209
 
          </relative_error>
210
 
          <max_iterations>
211
 
            <integer_value rank="0">500</integer_value>
212
 
          </max_iterations>
213
 
          <never_ignore_solver_failures/>
214
 
          <cache_solver_context/>
215
 
          <diagnostics>
216
 
            <monitors/>
217
 
          </diagnostics>
218
 
        </solver>
219
 
        <initial_condition name="WholeMesh">
220
 
          <python>
221
 
            <string_value lines="20" type="code" language="python">def val(X, t):
222
 
  import shallow_water_optimisation as constants
223
 
  #return constants.eta_exact(X[0], 0.0)
224
 
  # Instead of prescribing the exact eta we will optimise to it
225
 
  return 0.0</string_value>
226
 
          </python>
227
 
        </initial_condition>
228
 
        <mean_layer_thickness>
229
 
          <real_value rank="0">1</real_value>
230
 
        </mean_layer_thickness>
231
 
        <output/>
232
 
        <stat/>
233
 
        <consistent_interpolation/>
234
 
      </prognostic>
235
 
    </scalar_field>
236
 
    <scalar_field name="AnalyticalLayerThickness" rank="0">
237
 
      <prescribed>
238
 
        <mesh name="PressureMesh"/>
239
 
        <value name="WholeMesh">
240
 
          <python>
241
 
            <string_value lines="20" type="code" language="python">def val(X,t):
242
 
 import shallow_water_optimisation as constants
243
 
 x=X[0]
244
 
 return constants.eta_exact(x, t)</string_value>
245
 
          </python>
246
 
        </value>
247
 
        <output/>
248
 
        <stat/>
249
 
        <detectors>
250
 
          <exclude_from_detectors/>
251
 
        </detectors>
252
 
        <adjoint_storage>
253
 
          <exists_in_forward/>
254
 
        </adjoint_storage>
255
 
      </prescribed>
256
 
    </scalar_field>
257
 
    <scalar_field name="LayerThicknessError" rank="0">
258
 
      <diagnostic>
259
 
        <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">
260
 
          <absolute_difference/>
261
 
        </algorithm>
262
 
        <mesh name="PressureMesh"/>
263
 
        <output/>
264
 
        <stat/>
265
 
        <convergence>
266
 
          <include_in_convergence/>
267
 
        </convergence>
268
 
        <detectors>
269
 
          <include_in_detectors/>
270
 
        </detectors>
271
 
        <steady_state>
272
 
          <include_in_steady_state/>
273
 
        </steady_state>
274
 
        <adjoint_storage>
275
 
          <exists_in_forward/>
276
 
        </adjoint_storage>
277
 
      </diagnostic>
278
 
    </scalar_field>
279
 
    <vector_field name="AnalyticalVelocity" rank="1">
280
 
      <prescribed>
281
 
        <mesh name="VelocityMesh"/>
282
 
        <value name="WholeMesh">
283
 
          <python>
284
 
            <string_value lines="20" type="code" language="python">def val(X,t):
285
 
 import shallow_water_optimisation as constants
286
 
 x = X[0]
287
 
 return constants.u_exact(x, t)</string_value>
288
 
          </python>
289
 
        </value>
290
 
        <output/>
291
 
        <stat>
292
 
          <include_in_stat/>
293
 
        </stat>
294
 
        <detectors>
295
 
          <exclude_from_detectors/>
296
 
        </detectors>
297
 
        <adjoint_storage>
298
 
          <exists_in_forward/>
299
 
        </adjoint_storage>
300
 
      </prescribed>
301
 
    </vector_field>
302
 
    <vector_field name="VelocityError" rank="1">
303
 
      <diagnostic>
304
 
        <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">
305
 
          <absolute_difference/>
306
 
        </algorithm>
307
 
        <mesh name="VelocityMesh"/>
308
 
        <output/>
309
 
        <stat>
310
 
          <include_in_stat/>
311
 
        </stat>
312
 
        <convergence>
313
 
          <include_in_convergence/>
314
 
        </convergence>
315
 
        <detectors>
316
 
          <include_in_detectors/>
317
 
        </detectors>
318
 
        <steady_state>
319
 
          <include_in_steady_state/>
320
 
        </steady_state>
321
 
        <adjoint_storage>
322
 
          <exists_in_forward/>
323
 
        </adjoint_storage>
324
 
      </diagnostic>
325
 
    </vector_field>
326
 
  </material_phase>
327
 
  <adjoint>
328
 
    <functional name="integral_eta_t1">
329
 
      <functional_value>
330
 
        <algorithm name="functional_value">
331
 
          <string_value lines="20" type="code" language="python">J = 0.0
332
 
T1 = 0.5 # the time at which to evaluate
333
 
T2 = 1.0 # the time at which to evaluate
334
 
T = -1
335
 
if time &lt; T1 &lt;= time+dt:
336
 
  T = T1
337
 
elif time &lt; T2 &lt;= time+dt:
338
 
  T = T2
339
 
if time &lt; T &lt;= time+dt:
340
 
  import numpy
341
 
  eta_prev = states[n-1]["Fluid"].scalar_fields["LayerThickness"]
342
 
  eta      = states[n]["Fluid"].scalar_fields["LayerThickness"]
343
 
  eta_exact_prev = states[n-1]["Fluid"].scalar_fields["AnalyticalLayerThickness"]
344
 
  eta_exact      = states[n]["Fluid"].scalar_fields["AnalyticalLayerThickness"]
345
 
    
346
 
  # We want to temporally interpolate to evaluate eta at t=1.0
347
 
  alpha = (time + dt - T) / dt
348
 
  assert 0 &lt;= alpha &lt; 1
349
 
  tmp_eta = alpha * (eta_prev.val - eta_exact_prev.val) + (1-alpha) * (eta.val - eta_exact.val)
350
 
  
351
 
  # Now we want to integrate that over space
352
 
  coord = states[0]["Fluid"].vector_fields["Coordinate"]
353
 
  assert eta.element_count == eta_prev.element_count == coord.element_count
354
 
  for ele in range(coord.element_count):
355
 
    t = Transform(ele, coord)
356
 
    shape = eta_prev.ele_shape(ele)
357
 
    mass = t.shape_shape(shape, shape)
358
 
    nodes = eta_prev.ele_nodes(ele)
359
 
    J = J + numpy.dot(tmp_eta[nodes], numpy.dot(mass, tmp_eta[nodes]))</string_value>
360
 
        </algorithm>
361
 
        <reduction>
362
 
          <sum/>
363
 
        </reduction>
364
 
      </functional_value>
365
 
      <functional_dependencies>
366
 
        <algorithm name="functional_dependencies">
367
 
          <string_value lines="20" type="code" language="python">def dependencies(times, timestep):
368
 
  if times[0] &lt; 0.5 &lt;= times[1] or times[0] &lt; 1.0 &lt;= times[1]:
369
 
    return {"Fluid::Coordinate": [0],
370
 
            "Fluid::LayerThickness": [timestep-1, timestep],
371
 
            "Fluid::AnalyticalLayerThickness": [timestep-1, timestep]
372
 
            }
373
 
  else:
374
 
    return {}</string_value>
375
 
        </algorithm>
376
 
      </functional_dependencies>
377
 
    </functional>
378
 
    <controls>
379
 
      <control name="InitEta">
380
 
        <type field_name="Fluid::LayerThickness" name="initial_condition"/>
381
 
      </control>
382
 
      <load_controls/>
383
 
    </controls>
384
 
    <debug>
385
 
      <html_output/>
386
 
    </debug>
387
 
  </adjoint>
388
 
</shallow_water_options>