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