← Back to branch summary

~reducedmodelling/fluidity/ROM_Non-intrusive-ann

« back to all changes in this revision

Viewing changes to enkf/supermesh.F90-5Nov12

  • Committer: fangf at ac
  • Date: 2012-11-06 12:21:31 UTC
  • mto: This revision was merged to the branch mainline in revision 3989.
  • Revision ID: fangf@imperial.ac.uk-20121106122131-u2zvt7fxc1r3zeou
updated

Show diffs side-by-side

added added

removed removed

Lines of Context:
enkf/supermesh.F90-5Nov12
 
1
#include "fdebug.h"
 
2
 
 
3
module interpolation_ensemble_state_on_supermesh
 
4
  use fields
 
5
  use state_module
 
6
  use vtk_interfaces
 
7
  use pseudo_supermesh
 
8
  use interpolation_manager
 
9
  use spud
 
10
  use unify_meshes_module
 
11
  use mesh_files
 
12
  use populate_state_module
 
13
  use populate_sub_state_module
 
14
  use Profiler
 
15
 
 
16
 
 
17
 
 
18
   use elements
 
19
  use state_module
 
20
  use FLDebug
 
21
!  use spud
 
22
  use mesh_files
 
23
  use vtk_cache_module
 
24
  use global_parameters, only: OPTION_PATH_LEN, is_active_process, pi, &
 
25
    no_active_processes, topology_mesh_name, adaptivity_mesh_name, &
 
26
    periodic_boundary_option_path, domain_bbox, domain_volume
 
27
  use field_options
 
28
  use reserve_state_module
 
29
  use fields_manipulation
 
30
  use diagnostic_variables, only: convergence_field, steady_state_field
 
31
  use field_options
 
32
  use surfacelabels
 
33
  use climatology
 
34
  use metric_tools
 
35
  use coordinates
 
36
  use halos
 
37
  use tictoc
 
38
  use hadapt_extrude
 
39
  use hadapt_extrude_radially
 
40
  use initialise_fields_module
 
41
  use transform_elements
 
42
  use parallel_tools
 
43
  use boundary_conditions_from_options
 
44
  use nemo_states_module
 
45
  use data_structures
 
46
  use fields_halos
 
47
  use read_triangle
 
48
  use sediment, only: get_nSediments, get_sediment_name
 
49
 implicit none
 
50
 
 
51
contains
 
52
 
 
53
  subroutine interpolation_ensembles_on_supermesh(ensemble_state_new, ensemble_state_old)
 
54
 
 
55
  type(state_type), dimension(:), intent(inout) :: ensemble_state_old
 
56
  type(state_type), dimension(:), intent(inout) :: ensemble_state_new
 
57
 
 
58
 
 
59
  type(state_type) :: tmp_state_1, tmp_state_2
 
60
 
 
61
  type(vector_field), dimension(:), allocatable :: initial_positions
 
62
  type(vector_field) :: positions_old,positions_new,out_positions
 
63
  !type(vector_field),pointer :: out_positions
 
64
  character(len=255) :: filename
 
65
  character(len=255), dimension(:), allocatable :: files
 
66
  integer :: argc
 
67
  integer :: i, status,ifield,nfield
 
68
  type(state_type) :: initial_state
 
69
  integer :: ierr
 
70
  integer :: stat, mxnods
 
71
  integer :: nrens
 
72
  integer :: field_count
 
73
  type(vector_field) :: velocity
 
74
  type(scalar_field) :: pressure, free_surface,time
 
75
  type(element_type) :: shape
 
76
  integer :: quad_degree
 
77
 
 
78
  type(vector_field),pointer :: field
 
79
   type(vector_field):: vfield2
 
80
   type(scalar_field):: sfield2
 
81
   type(tensor_field):: tfield2
 
82
 
 
83
   type(vector_field), pointer :: vfield,vfield_old,vfield_new
 
84
   type(scalar_field), pointer :: sfield,sfield_old,sfield_new
 
85
   type(tensor_field), pointer :: tfield,tfield_old,tfield_new
 
86
 
 
87
 
 
88
    print*,'In interpolation_ensembles_on_supermesh'
 
89
    nrens = size(ensemble_state_old)    
 
90
    mxnods = 100000
 
91
 
 
92
    print*,nrens    
 
93
 
 
94
    call mpi_init(ierr)
 
95
    call set_option('/mesh_adaptivity/hr_adaptivity/maximum_number_of_nodes', mxnods, stat=stat)
 
96
    call get_option('/geometry/quadrature/degree', quad_degree)    
 
97
 
 
98
 
 
99
!+++++++++++++++++++++++++++++++++++++
 
100
    allocate(initial_positions(nrens))
 
101
    do i=1,nrens
 
102
       initial_positions(i) = extract_vector_field(ensemble_state_old(i), "Coordinate")
 
103
    enddo
 
104
 
 
105
    allocate(files(nrens))
 
106
    !do i=1, nrens
 
107
    !   !! Note that this won't work in parallel. Have to look for the pvtu in that case.
 
108
    !   write(files(i), '(a, i0, a)') trim(simulation_name)//'_', i, ".vtu"
 
109
    !end do
 
110
    
 
111
    !goto 300
 
112
    !call vtk_read_state(trim(files(1)), initial_state)
 
113
    !initial_positions = extract_vector_field(initial_state, "Coordinate")
 
114
    !call add_faces(initial_positions%mesh)
 
115
    
 
116
    !goto 200 
 
117
    print*,trim(files(1)),trim(files(2))
 
118
    !!Generate supermesh
 
119
    !call compute_pseudo_supermesh(files, initial_positions(1), out_positions, mxnods=mxnods)
 
120
    
 
121
    !initial_positions(1)=read_mesh_files(files(1), 'gmsh', shape)
 
122
    !initial_positions(1)=read_mesh_files('munk_gyre-ns_checkpoint_checkpoint_CoordinateMesh_1', 'gmsh', shape)
 
123
    initial_positions(1)=read_mesh_files('munk_gyre-ns_checkpoint_checkpoint_CoordinateMesh_1', quad_degree=quad_degree)
 
124
    initial_positions(2)=read_mesh_files('munk_gyre-ns_checkpoint_checkpoint_CoordinateMesh_2', quad_degree=quad_degree)
 
125
    
 
126
    
 
127
    print*,trim(initial_positions(1)%mesh%name)
 
128
    
 
129
    print*,'before compute_pseudo_supermesh'
 
130
    call compute_pseudo_supermesh_from_mesh(files, initial_positions(:), out_positions, mxnods=mxnods)
 
131
    
 
132
    
 
133
    do i=1, nrens
 
134
       
 
135
       
 
136
       call insert(ensemble_state_new(i), out_positions, "Coordinate")
 
137
       call insert(ensemble_state_new(i), out_positions%mesh, "CoordinateMesh")
 
138
    end do
 
139
    
 
140
    print*,trim(initial_positions(1)%mesh%name)
 
141
    
 
142
    !print*,'before compute_pseudo_supermesh'
 
143
    !call compute_pseudo_supermesh_from_mesh(files, initial_positions(:), out_positions, mxnods=mxnods)
 
144
    
 
145
    call insert_derived_meshes( ensemble_state_new)
 
146
    call allocate_and_insert_fields(ensemble_state_new)
 
147
    
 
148
    
 
149
    do i = 2, nrens
 
150
       nfield = scalar_field_count( ensemble_state_new(1) )
 
151
       do ifield = 1, nfield 
 
152
          sfield => extract_scalar_field( ensemble_state_new(1), ifield )
 
153
          call allocate( sfield2,  sfield%mesh, trim(sfield%name)   )
 
154
          call set(sfield2, sfield)
 
155
          call insert(ensemble_state_new(i), sfield2, trim(sfield%name))
 
156
       end do
 
157
       
 
158
       
 
159
       nfield = vector_field_count( ensemble_state_new(1) )
 
160
       print*,'&&&&vector',nfield
 
161
       do ifield = 1, nfield 
 
162
          vfield => extract_vector_field( ensemble_state_new(1), ifield )
 
163
          print *, node_count(vfield),trim(vfield%name)
 
164
          call allocate( vfield2,  vfield%dim, vfield%mesh, trim(vfield%name)   )
 
165
          call set(vfield2, vfield)
 
166
          call insert(ensemble_state_new(i), vfield2, trim(vfield%name))
 
167
       end do
 
168
       
 
169
       nfield = tensor_field_count( ensemble_state_new(1) )
 
170
       do ifield = 1, nfield 
 
171
          tfield => extract_tensor_field( ensemble_state_new(1), ifield )
 
172
          call allocate( tfield2,  tfield%mesh, trim(tfield%name)    )
 
173
          call set(tfield2, tfield)
 
174
          call insert(ensemble_state_new(i), tfield2, trim(tfield%name))
 
175
       end do
 
176
    end do
 
177
    
 
178
    
 
179
    field=> extract_vector_field(ensemble_state_new(1), "Velocity")
 
180
    print *, node_count(field)
 
181
    !    print*,'before unify_meshes'
 
182
    !    out_positions=unify_meshes(initial_positions)
 
183
    
 
184
    print*,node_count(initial_positions(1)%mesh)
 
185
    print*,node_count(initial_positions(2)%mesh)
 
186
    print*,node_count(out_positions%mesh)
 
187
    
 
188
    
 
189
    !! call vtk_write_fields("pseudo_supermesh", 0, out_positions, out_positions%mesh)
 
190
    
 
191
    !! Double check whether ensemble_state_new is setup correctly?? miss fields???
 
192
    print*,'after allocate pressure'
 
193
    do i=1, nrens
 
194
       call print_state( ensemble_state_new(i) )
 
195
       print *, 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
 
196
       call print_state( ensemble_state_old(i) )
 
197
       
 
198
    end do
 
199
    
 
200
    !+++++++++++++++++++++++++++++++++++++
 
201
    
 
202
    print*,'linear_interpolation'
 
203
    do i=1, nrens
 
204
       print*,'****iren=',i
 
205
       print*,'interpolation of scalar'
 
206
       nfield = scalar_field_count( ensemble_state_old(i) )
 
207
       do ifield = 1, nfield 
 
208
          sfield_old => extract_scalar_field( ensemble_state_old(i), ifield )
 
209
          !    call allocate( sfield_old,  sfield%mesh, trim(sfield%name)   )
 
210
          print*,trim(sfield_old%name)
 
211
          !    call set(sfield_old, sfield)
 
212
          positions_old = extract_vector_field(ensemble_state_old(i), "Coordinate")
 
213
          ! new field
 
214
          sfield_new => extract_scalar_field( ensemble_state_new(i), ifield )
 
215
          !       call allocate( sfield_new,  sfield%mesh, trim(sfield%name)   )
 
216
          print*,trim(sfield_new%name)
 
217
          !       call set(sfield_new, sfield)
 
218
          positions_new = extract_vector_field(ensemble_state_new(i), "Coordinate")
 
219
          call linear_interpolation(sfield_old, positions_old, sfield_new,positions_new)
 
220
       end do
 
221
       
 
222
       
 
223
       print*,'interpolation of vector'
 
224
       nfield = vector_field_count( ensemble_state_new(i) )
 
225
       do ifield = 1, nfield-1 
 
226
          print*,'ifield',ifield
 
227
          vfield_old => extract_vector_field( ensemble_state_old(i), ifield )
 
228
          !          call allocate( vfield_old, vfield%dim, vfield%mesh, trim(vfield%name)   )
 
229
          print*,trim(vfield_old%name)
 
230
          print*,'11'
 
231
          !          call set(vfield_old, vfield)
 
232
          positions_old = extract_vector_field(ensemble_state_old(i), "Coordinate")
 
233
          ! new field
 
234
          vfield_new => extract_vector_field( ensemble_state_new(i), ifield )
 
235
          !  call allocate( vfield_new, vfield%dim,vfield%mesh, trim(vfield%name)   )
 
236
          print*,trim(vfield_new%name)
 
237
          print*,'22'
 
238
          !  call set(vfield_new, vfield)
 
239
          positions_new = extract_vector_field(ensemble_state_new(i), "Coordinate")
 
240
          call linear_interpolation(vfield_old, positions_old, vfield_new,positions_new)
 
241
       end do
 
242
       
 
243
       print*,'interpolation of tensor'
 
244
       nfield = tensor_field_count( ensemble_state_new(i) )
 
245
       do ifield = 1, nfield 
 
246
          tfield_old => extract_tensor_field( ensemble_state_old(i), ifield )
 
247
          !    call allocate( tfield_old,  tfield%mesh, trim(tfield%name)   )
 
248
          !    call set(tfield_old, tfield)
 
249
          positions_old = extract_vector_field(ensemble_state_old(i), "Coordinate")
 
250
          ! new field
 
251
          tfield_new => extract_tensor_field( ensemble_state_new(i), ifield )
 
252
          !   call allocate( tfield_new,  tfield%mesh, trim(tfield%name)   )
 
253
          !   call set(tfield_new, tfield)
 
254
          positions_new = extract_vector_field(ensemble_state_new(i), "Coordinate")
 
255
          call linear_interpolation(tfield_old, positions_old, tfield_new,positions_new)
 
256
          
 
257
       end do
 
258
       print*,'%%%%%%%%%%%%%%%%%%%%%%%%%i=',i
 
259
       !       call linear_interpolation(ensemble_state_old(i), ensemble_state_new(i))
 
260
       pressure = extract_scalar_field(ensemble_state_new(i), "Pressure")
 
261
       print*,pressure%val(1:20)
 
262
       print*,'pressurenew',node_count(pressure%mesh)
 
263
       !      stop 23
 
264
    enddo
 
265
    
 
266
 
 
267
    print *, '55'
 
268
    stop 666
 
269
 
 
270
  end subroutine interpolation_ensembles_on_supermesh
 
271
 
 
272
 
 
273
    !!femtools/Fields_Calculations.F90
 
274
    !function norm2_difference_single(fieldA, positionsA, fieldB, positionsB) result(norm)
 
275
 
 
276
 
 
277
  subroutine compute_pseudo_supermesh_from_mesh(snapshots, starting_positions, super_positions, no_its, mxnods)
 
278
    !!< snapshots is a list of VTUs containing the meshes we want
 
279
    !!< to merge.
 
280
    !!< starting_positions is the initial mesh + positions to interpolate the
 
281
    !!< metric tensor fields describing the snapshot meshes onto.
 
282
    !!< super_positions is the output -- a positions field on a mesh.
 
283
    character(len=255), dimension(:), intent(in) :: snapshots
 
284
    type(vector_field), dimension(:), intent(in) :: starting_positions
 
285
    type(vector_field), intent(out) :: super_positions
 
286
    integer, intent(in), optional :: no_its, mxnods
 
287
 
 
288
    integer :: lno_its
 
289
    integer :: it, i
 
290
 
 
291
    type(mesh_type) :: current_mesh, vtk_mesh
 
292
    type(vector_field) :: current_pos, vtk_pos
 
293
    type(state_type) :: vtk_state, temp_state
 
294
    type(state_type) :: interpolation_input, interpolation_output
 
295
 
 
296
    type(tensor_field) :: merged_metric, interpolated_metric
 
297
    type(tensor_field) :: vtk_metric
 
298
 
 
299
 
 
300
    if (present(no_its)) then
 
301
      lno_its = no_its
 
302
    else
 
303
      lno_its = 3
 
304
    end if
 
305
 
 
306
    current_pos = starting_positions(1)
 
307
    call incref(starting_positions(1))
 
308
    current_mesh = starting_positions(1)%mesh
 
309
    call incref(current_mesh)
 
310
 
 
311
    do it=1,lno_its
 
312
      call allocate(merged_metric, current_mesh, "MergedMetric")
 
313
      call zero(merged_metric)
 
314
      call allocate(interpolated_metric, current_mesh, "InterpolatedMetric")
 
315
      call zero(interpolated_metric)
 
316
      call insert(interpolation_output, interpolated_metric, "InterpolatedMetric")
 
317
      call insert(interpolation_output, current_mesh, "CoordinateMesh")
 
318
      call insert(interpolation_output, current_pos, "Coordinate")
 
319
 
 
320
!      call allocate(edgelen, current_mesh, "EdgeLengths")
 
321
 
 
322
      do i=1,size(snapshots)
 
323
        call zero(interpolated_metric)
 
324
        !!call vtk_read_state(trim(snapshots(i)), vtk_state)
 
325
        !!vtk_mesh = extract_mesh(vtk_state, "Mesh")
 
326
        !!vtk_pos  = extract_vector_field(vtk_state, "Coordinate")
 
327
        vtk_mesh = starting_positions(i)%mesh
 
328
        vtk_pos  = starting_positions(i)
 
329
 
 
330
        call allocate(vtk_metric, vtk_mesh, "MeshMetric")
 
331
        print*,'before compute_mesh_metric'
 
332
        call compute_mesh_metric(vtk_pos, vtk_metric)
 
333
        print*,'after compute_mesh_metric'
 
334
        call insert(interpolation_input, vtk_metric, "InterpolatedMetric")
 
335
        call insert(interpolation_input, vtk_mesh, "CoordinateMesh")
 
336
        call insert(interpolation_input, vtk_pos, "Coordinate")
 
337
        print*,'vtk_write_state("interpolation_input")'
 
338
        !!call vtk_write_state("interpolation_input", i, state=(/interpolation_input/))
 
339
        call linear_interpolation(interpolation_input, interpolation_output)
 
340
        print*,'vtk_write_state("interpolation_output")'
 
341
        !!call vtk_write_state("interpolation_output", i, state=(/interpolation_output/))
 
342
        call merge_tensor_fields(merged_metric, interpolated_metric)
 
343
        call deallocate(vtk_metric)
 
344
        call deallocate(interpolation_input)
 
345
        call deallocate(vtk_state)
 
346
      end do
 
347
 
 
348
      call deallocate(interpolated_metric)
 
349
      call deallocate(interpolation_output)
 
350
 
 
351
      call insert(temp_state, current_mesh, "CoordinateMesh")
 
352
      call insert(temp_state, current_pos, "Coordinate")
 
353
      ! Assuming current_mesh had a refcount of one,
 
354
      ! it now has a refcount of two.
 
355
!      call get_edge_lengths(merged_metric, edgelen)
 
356
!      call vtk_write_fields("supermesh_before_adapt", it, current_pos, current_mesh, sfields=(/edgelen/), tfields=(/merged_metric/))
 
357
      if (present(mxnods)) then
 
358
        call limit_metric(current_pos, merged_metric, min_nodes=1, max_nodes=mxnods)
 
359
      end if
 
360
      print*,'before adapt_state'
 
361
      call adapt_state(temp_state, merged_metric)
 
362
!      call vtk_write_state("supermesh_after_adapt", it, state=(/temp_state/))
 
363
      call deallocate(merged_metric)
 
364
!      call deallocate(edgelen)
 
365
      ! Now it has a refcount of one, as adapt_state
 
366
      ! has destroyed the old one and created a new mesh
 
367
      ! with refcount one.
 
368
 
 
369
      ! We're finished with the current_mesh, so let it be
 
370
      ! deallocated if no one else is using it.
 
371
      call deallocate(current_mesh)
 
372
      call deallocate(current_pos)
 
373
 
 
374
      current_mesh = extract_mesh(temp_state, "CoordinateMesh")
 
375
      current_pos  = extract_vector_field(temp_state, "Coordinate")
 
376
      call incref(current_mesh)
 
377
      call incref(current_pos)
 
378
      call deallocate(temp_state)
 
379
    end do
 
380
 
 
381
    call deallocate(current_mesh)
 
382
    super_positions = current_pos
 
383
  end subroutine compute_pseudo_supermesh_from_mesh
 
384
 
 
385
end module Interpolation_ensemble_state_on_supermesh