SFEMaNS  version 5.3
Reference documentation for SFEMaNS
post_processing_debug.f90
Go to the documentation of this file.
1 MODULE post_processing_debug
2 
3 PUBLIC :: post_proc_test
4 
5 PRIVATE
6 CONTAINS
7  !---------------------------------------------------------------------------
8  SUBROUTINE post_proc_test(vv_mesh, pp_mesh, temp_mesh, H_mesh, phi_mesh, list_mode, &
9  un, pn, hn, bn, phin, tempn, level_setn, mu_h_field, &
10  time, m_max_c, comm_one_d, comm_one_d_ns, comm_one_d_temp) ! MODIFICATION: comm_one_d_temp added
11  USE boundary
12  USE def_type_mesh
13  USE input_data
14  USE my_util
15  USE tn_axi
16  USE subroutine_ns_with_u
17  USE sft_parallele
18 #include "petsc/finclude/petsc.h"
19  USE petsc
20  IMPLICIT NONE
21  TYPE(mesh_type), POINTER :: pp_mesh, vv_mesh
22  TYPE(mesh_type), POINTER :: temp_mesh
23  TYPE(mesh_type), POINTER :: H_mesh, phi_mesh
24  INTEGER, POINTER, DIMENSION(:) :: list_mode
25  REAL(KIND=8), POINTER, DIMENSION(:,:,:) :: un, pn, Hn, Bn, phin, tempn
26  REAL(KIND=8), POINTER, DIMENSION(:,:,:,:) :: level_setn
27  REAL(KIND=8), POINTER, DIMENSION(:) :: mu_H_field
28  REAL(KIND=8) :: time
29  INTEGER :: m_max_c
30  REAL(KIND=8), DIMENSION(SIZE(un,1), SIZE(un,2), SIZE(un,3)) :: un_m1, un_ex, un_error
31  REAL(KIND=8), DIMENSION(SIZE(pn,1), SIZE(pn,2), SIZE(pn,3)) :: pn_m1, pn_ex, pn_error
32  REAL(KIND=8), DIMENSION(SIZE(Hn,1), SIZE(Hn,2), SIZE(Hn,3)) :: Hn1, Hn_ex, Hn_error
33  REAL(KIND=8), DIMENSION(SIZE(phin,1),SIZE(phin,2),SIZE(phin,3)) :: phin1
34  REAL(KIND=8), DIMENSION(SIZE(tempn,1), SIZE(tempn,2), SIZE(tempn,3)) :: tempn_m1, tempn_ex, tempn_error
35  REAL(KIND=8), DIMENSION(SIZE(level_setn,1),SIZE(level_setn,2),SIZE(level_setn,3),SIZE(level_setn,4)):: level_setn_m1
36 
37  INTEGER :: i, k, code, rank, int_nb, n, TYPE
38  REAL(KIND=8), DIMENSION(4) :: norm_err_loc, norm_err
39  REAL(KIND=8) :: err, norm
40  REAL(KIND=8) :: moyenne
41  mpi_comm, DIMENSION(:), POINTER :: comm_one_d, comm_one_d_ns, comm_one_d_temp ! MODIFICATION: comm_one_d_temp added
42 
43  CALL mpi_comm_rank(mpi_comm_world,rank,code)
44 
45  SELECT CASE(inputs%numero_du_test_debug)
46  CASE(1,2,8,9,19,20,24,25,39,40)
47  DO i = 1, m_max_c
48  DO k= 1, 6
49  un_m1(:,k,i) = un(:,k,i) - vv_exact(k,vv_mesh%rr,list_mode(i),time)
50  END DO
51  DO k= 1, 2
52  pn_m1(:,k,i) = pn(:,k,i) - pp_exact(k,pp_mesh%rr,list_mode(i),time)
53  END DO
54  IF (list_mode(i) == 0) THEN
55  CALL moy(comm_one_d(1),pp_mesh, pn_m1(:,1,i),moyenne)
56  pn_m1(:,1,i) = pn_m1(:,1,i) - moyenne
57  ENDIF
58  END DO
59 
60  !norm_err(1) = norm_SF(comm_one_d_NS, 'L2', vv_mesh, list_mode, un_m1)
61  norm_err(1) = sqrt(dot_product_sf(comm_one_d_ns,vv_mesh, list_mode, un_m1, un_m1))
62  norm_err(2) = norm_sf(comm_one_d_ns, 'sH1', vv_mesh, list_mode, un_m1)
63  norm_err(3) = norm_sf(comm_one_d_ns, 'div', vv_mesh, list_mode, un)
64  norm_err(4) = norm_sf(comm_one_d_ns, 'L2', pp_mesh, list_mode, pn_m1)
65  IF (rank==0) THEN
66  WRITE(10,*) 'Velocity field #####################'
67  WRITE(10,*) 'L2 error on velocity = ', norm_err(1)
68  WRITE(10,*) 'H1 error on velocity = ', norm_err(2)
69  WRITE(10,*) 'L2 norm of divergence = ', norm_err(3)
70  WRITE(10,*) 'Pressure field #####################'
71  WRITE(10,*) 'L2 error on pressure = ', norm_err(4)
72  END IF
73 
74  IF (inputs%if_temperature) THEN
75  DO i = 1, m_max_c
76  DO k= 1, 2
77  tempn_m1(:,k,i) = tempn(:,k,i) - temperature_exact(k,temp_mesh%rr, list_mode(i), time)
78  END DO
79  END DO
80  norm_err(1) = norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_m1) ! MODIFICATION: comm_one_d_temp instead of ns
81  norm_err(2) = norm_sf(comm_one_d_temp, 'sH1', temp_mesh, list_mode, tempn_m1) ! MODIFICATION: comm_one_d_temp instead of ns
82  IF (rank==0) THEN
83  WRITE(10,*) 'Temperature field#####################'
84  WRITE(10,*) 'L2 error on temperature = ', norm_err(1)
85  WRITE(10,*) 'H1 error on temperature = ', norm_err(2)
86  WRITE(10,*) 'Temperature field#####################'
87  END IF
88  END IF
89 
90  IF (inputs%if_level_set) THEN
91  IF (inputs%if_level_set_P2) THEN
92  DO i = 1, m_max_c
93  DO k = 1, 2
94  DO int_nb = 1, inputs%nb_fluid - 1
95  level_setn_m1(int_nb,:,k,i) = level_setn(int_nb,:,k,i) &
96  -level_set_exact(int_nb,k,vv_mesh%rr,list_mode(i),time)
97  END DO
98  END DO
99  END DO
100  IF (inputs%numero_du_test_debug==39) THEN
101  norm_err(3) = norm_s_l1_zero_mode(comm_one_d_ns, vv_mesh, list_mode,level_setn_m1(1,:,:,:))
102  IF (rank==0) THEN
103  WRITE(10,*) 'Level set field#####################'
104  WRITE(10,*) 'L1 error on level set', norm_err(3)
105  WRITE(10,*) 'Level set field#####################'
106  END IF
107  ELSE
108  norm_err(3) = norm_sf(comm_one_d_ns, 'L2', vv_mesh, list_mode,level_setn_m1(1,:,:,:))
109  IF (rank==0) THEN
110  WRITE(10,*) 'Level set field#####################'
111  WRITE(10,*) 'L2 error on level set', norm_err(3)
112  WRITE(10,*) 'Level set field#####################'
113  END IF
114  END IF
115  ELSE
116  DO i = 1, m_max_c
117  DO k = 1, 2
118  DO int_nb = 1, inputs%nb_fluid - 1
119  level_setn_m1(int_nb,:,k,i) = level_setn(int_nb,:,k,i) &
120  -level_set_exact(int_nb,k,pp_mesh%rr,list_mode(i),time)
121  END DO
122  END DO
123  END DO
124  IF (inputs%numero_du_test_debug==39) THEN
125  norm_err(3) = norm_s_l1_zero_mode(comm_one_d_ns, pp_mesh, list_mode,level_setn_m1(1,:,:,:))
126  IF (rank==0) THEN
127  WRITE(10,*) 'Level set field#####################'
128  WRITE(10,*) 'L1 error on level set', norm_err(3)
129  WRITE(10,*) 'Level set field#####################'
130  END IF
131  ELSE
132  norm_err(3) = norm_sf(comm_one_d_ns, 'L2', pp_mesh, list_mode,level_setn_m1(1,:,:,:))
133  IF (rank==0) THEN
134  WRITE(10,*) 'Level set field#####################'
135  WRITE(10,*) 'L2 error on level set', norm_err(3)
136  WRITE(10,*) 'Level set field#####################'
137  END IF
138  END IF
139  END IF
140  END IF
141 
142  IF (inputs%numero_du_test_debug==24) THEN
143  DO i = 1, m_max_c
144  DO k= 1, 2
145  DO n = 1, SIZE(pp_mesh%rr,2)
146  IF (pp_mesh%rr(1,n).GT..5d0) THEN
147  pn_m1(n,k,i) = pn_m1(n,k,i)
148  ELSE
149  pn_m1(n,k,i) = 0.d0
150  END IF
151  END DO
152  END DO
153  END DO
154 
155  norm_err(4) = norm_sf(comm_one_d_ns, 'L2', pp_mesh, list_mode, pn_m1)
156  IF (rank==0) THEN
157  WRITE(10,*) 'L2 error on pressure outter obstacle = ', norm_err(4)
158  END IF
159  END IF
160 
161 
162  CASE(3:7,10,12)
163  DO i = 1, m_max_c
164  DO k =1, 6
165  hn1(:,k,i) = hn(:,k,i) - hexact(h_mesh, k, h_mesh%rr, list_mode(i), mu_h_field, time)
166  END DO
167  END DO
168 
169  IF (inputs%nb_dom_phi/=0) THEN
170  DO i = 1, m_max_c
171  DO k =1, 2
172  phin1(:,k,i) = phin(:,k,i) - phiexact(k, phi_mesh%rr, list_mode(i), inputs%mu_phi, time)
173  END DO
174  END DO
175  END IF
176 
177  norm = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn)
178  err = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn1)
179  norm_err(1) = err/norm
180 
181  norm = norm_sf(comm_one_d, 'sH1', h_mesh, list_mode, hn)
182  err = norm_sf(comm_one_d, 'curl', h_mesh, list_mode, hn1)
183  norm_err(2) = err/norm
184 
185  norm = norm_sf(comm_one_d, 'sH1', h_mesh, list_mode, bn)
186  err = norm_sf(comm_one_d, 'div', h_mesh, list_mode, bn)
187 
188  norm_err(3) = err/norm
189 
190  IF (inputs%nb_dom_phi/=0) THEN
191  norm = norm_sf(comm_one_d, 'sH1', phi_mesh, list_mode, phin)
192  err = norm_sf(comm_one_d, 'sH1', phi_mesh, list_mode, phin1)
193  norm_err(4) = err/norm
194  ELSE
195  norm_err(4) = inputs%norm_ref(4)
196  END IF
197 
198  IF (rank==0) THEN
199  WRITE(rank+10,*) 'Magnetic field #####################'
200  WRITE(rank+10,*) 'L2 norm of error on Hn = ', norm_err(1)
201  WRITE(rank+10,*) 'L2 norm of error on Curl(Hn) = ', norm_err(2)
202  WRITE(rank+10,*) 'L2 norm of Div(mu Hn) = ', norm_err(3)
203  IF (inputs%nb_dom_phi/=0) THEN
204  WRITE(10,*) 'Scal potential #####################'
205  WRITE(10,*) 'H1 norm of error on phin = ', norm_err(4)
206  END IF
207  END IF
208 
209  CASE(11)
210  norm_err(1) = norm_sf(comm_one_d_ns, 'sH1', vv_mesh, list_mode, un)
211  norm_err(2) = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn)
212  norm_err(3) = norm_sf(comm_one_d_ns, 'L2', pp_mesh, list_mode, pn)
213  norm_err(4) = norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn) ! MODIFICATION: comm_one_d_temp instead of ns
214  IF (rank==0) THEN
215  WRITE(10,*) '######################################'
216  WRITE(10,*) 'H1 norm on velocity = ', norm_err(1)
217  WRITE(10,*) 'L2 norm of Hn = ', norm_err(2)
218  WRITE(10,*) 'L2 norm of pressure = ', norm_err(3)
219  WRITE(10,*) 'L2 norm of temperature = ', norm_err(4)
220  WRITE(10,*) '######################################'
221  END IF
222 
223  CASE(13)
224  norm_err(1) = norm_sf(comm_one_d_ns, 'sH1', vv_mesh, list_mode, un)
225  norm_err(2) = norm_sf(comm_one_d, 'div', h_mesh, list_mode, hn)
226  norm_err(3) = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn)
227  norm_err(4) = norm_sf(comm_one_d_ns, 'L2', pp_mesh, list_mode, pn)
228  IF (rank==0) THEN
229  WRITE(10,*) '######################################'
230  WRITE(10,*) 'H1 norm on velocity = ', norm_err(1)
231  WRITE(10,*) 'L2 norm of div(Hn) = ', norm_err(2)
232  WRITE(10,*) 'L2 norm of Hn = ', norm_err(3)
233  WRITE(10,*) 'L2 norm of pressure = ', norm_err(4)
234  WRITE(10,*) '######################################'
235  END IF
236 
237  CASE(14)
238  IF (rank==0) THEN
239  READ(10,*) norm_err(1)
240  READ(10,*) norm_err(2)
241  READ(11,*) norm_err(3)
242  READ(11,*) norm_err(4)
243  END IF
244 
245  CASE(15)
246 
247  norm_err_loc = 0.d0
248  DO i = 1, m_max_c
249  norm_err_loc(rank+1)= 0.5*(norme_l2_champ_par(comm_one_d_ns(1), &
250  vv_mesh, list_mode(i:i), un(:,:,i:i)))**2
251  END DO
252  norm_err_loc(4)= norm_sf(comm_one_d_ns, 'div', vv_mesh, list_mode, un)/&
253  norm_sf(comm_one_d_ns, 'sH1', vv_mesh, list_mode, un)
254  CALL mpi_allreduce(norm_err_loc,norm_err,4,mpi_double_precision, &
255  mpi_max, comm_one_d(2), code)
256 
257  IF (rank==0) THEN
258  WRITE(10,*) '######################################'
259  WRITE(10,*) 'L2 norm of mode 0 = ', norm_err(1)
260  WRITE(10,*) 'L2 norm of mode 1 = ', norm_err(2)
261  WRITE(10,*) 'L2 norm of mode 2 = ', norm_err(3)
262  WRITE(10,*) 'Relative L2 norm of div.= ', norm_err(4)
263  WRITE(10,*) '######################################'
264  END IF
265 
266  CASE(16)
267  norm_err(1) = 0.5*norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un)**2
268  CALL angular_momentum(vv_mesh, list_mode, un, norm_err(2:4))
269  IF (rank==0) THEN
270  WRITE(10,*) '######################################'
271  WRITE(10,*) 'Total kinetic energy = ', norm_err(1)
272  WRITE(10,*) 'Mx = ', norm_err(2)
273  WRITE(10,*) 'My = ', norm_err(3)
274  WRITE(10,*) 'Mz = ', norm_err(4)
275  WRITE(10,*) '######################################'
276  END IF
277 
278  CASE(17,22,27)
279  norm = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, bn)
280  err = norm_sf(comm_one_d, 'div', h_mesh, list_mode, bn)
281  norm_err(1) = err/norm
282  DO i = 1, m_max_c
283  DO k = 1, 6
284  hn1(:,k,i) = hexact(h_mesh, k, h_mesh%rr, list_mode(i), mu_h_field, time)
285  hn1(:,k,i) = hn1(:,k,i) - hn(:,k,i)
286  ENDDO
287  ENDDO
288  err = norm_sf(comm_one_d, 'curl', h_mesh, list_mode, hn1)
289  norm = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn)
290  norm_err(2) = err/norm
291  err = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn1)
292  norm_err(3) = err/norm
293  DO i = 1, m_max_c
294  DO k = 1, 2
295  phin1(:,k,i) = phin(:,k,i) - phiexact(k, phi_mesh%rr, list_mode(i), inputs%mu_phi, time)
296  END DO
297  END DO
298  err = norm_sf(comm_one_d, 'L2', phi_mesh, list_mode, phin1)
299  norm = norm_sf(comm_one_d, 'L2', phi_mesh, list_mode, phin)
300  norm_err(4) = err/norm
301  IF (rank==0) THEN
302  WRITE(10,*) '########################################################'
303  WRITE(10,*) ' L2-norm on div of B err/norm =', norm_err(1)
304  WRITE(10,*) ' L2-norm on curl of (H-Hexact) err/norm =', norm_err(2)
305  WRITE(10,*) ' L2-norm of (H-Hexact) err/norm =', norm_err(3)
306  WRITE(10,*) ' L2-norm on phi err/norm =', norm_err(4)
307  WRITE(10,*) '########################################################'
308  END IF
309 
310  CASE(18,23,26)
311  norm = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, bn)
312  err = norm_sf(comm_one_d, 'div', h_mesh, list_mode, bn)
313  norm_err(1) = err/norm
314  DO i = 1, m_max_c
315  DO k = 1, 6
316  hn1(:,k,i) = hexact(h_mesh, k, h_mesh%rr, list_mode(i), mu_h_field, time)
317  hn1(:,k,i) = hn1(:,k,i) - hn(:,k,i)
318  ENDDO
319  ENDDO
320  err = norm_sf(comm_one_d, 'curl', h_mesh, list_mode, hn1)
321  norm = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn)
322  norm_err(2) = err/norm
323  err = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn1)
324  norm_err(3) = err/norm
325  norm_err(4) = norm
326  IF (rank==0) THEN
327  WRITE(10,*) '########################################################'
328  WRITE(10,*) ' L2-norm on div of B err/norm =', norm_err(1)
329  WRITE(10,*) ' L2-norm on curl of (H-Hexact) err/norm =', norm_err(2)
330  WRITE(10,*) ' L2-norm of (H-Hexact) err/norm =', norm_err(3)
331  WRITE(10,*) ' L2-norm of H =', norm_err(4)
332  WRITE(10,*) '########################################################'
333  END IF
334 
335  CASE(21)
336  DO i = 1, m_max_c
337  DO k= 1, 6
338  un_m1(:,k,i) = un(:,k,i) - vv_exact(k,vv_mesh%rr,list_mode(i),time)
339  hn1(:,k,i) = hn(:,k,i) - hexact(h_mesh, k, h_mesh%rr, list_mode(i), mu_h_field, time)
340  END DO
341  DO k= 1, 2
342  pn_m1(:,k,i) = pn(:,k,i) - pp_exact(k,pp_mesh%rr,list_mode(i),time)
343  DO int_nb = 1, inputs%nb_fluid - 1
344  level_setn_m1(int_nb,:,k,i) = level_setn(int_nb,:,k,i) &
345  - level_set_exact(int_nb,k,pp_mesh%rr,list_mode(i),time)
346  END DO
347  END DO
348  IF (list_mode(i) == 0) THEN
349  CALL moy(comm_one_d(1),pp_mesh, pn_m1(:,1,i),moyenne)
350  pn_m1(:,1,i) = pn_m1(:,1,i) - moyenne
351  ENDIF
352  END DO
353  norm_err(1) = norm_sf(comm_one_d_ns, 'L2', vv_mesh, list_mode, un_m1)
354  norm_err(2) = norm_sf(comm_one_d_ns, 'L2', pp_mesh, list_mode, pn_m1)
355  norm_err(3) = norm_sf(comm_one_d_ns, 'L2', pp_mesh, list_mode, level_setn_m1(1,:,:,:))
356  err = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn1)
357  norm_err(4) = err
358  IF (rank==0) THEN
359  WRITE(10,*) 'Velocity field #####################'
360  WRITE(10,*) 'L2 error on velocity = ', norm_err(1)
361  WRITE(10,*) 'Pressure field #####################'
362  WRITE(10,*) 'L2 error on pressure = ', norm_err(2)
363  WRITE(10,*) 'Level Set #####################'
364  WRITE(10,*) 'L2 error on level set = ', norm_err(3)
365  WRITE(10,*) 'Magnetic field #####################'
366  WRITE(10,*) 'L2 error on Hn = ', norm_err(4)
367  END IF
368 
369 
370  CASE(28)
371  err = norm_sf(comm_one_d, 'div', vv_mesh, list_mode, un)
372  norm = norm_sf(comm_one_d, 'H1', vv_mesh, list_mode, un)
373  norm_err(1) = err
374  norm_err(2) = err/norm
375  norm = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un)
376  norm_err(3) = norm
377  norm = norm_sf(comm_one_d, 'H1', pp_mesh, list_mode, pn)
378  norm_err(4) = norm
379 
380  IF (rank==0) THEN
381  WRITE(10,*) '########################################################'
382  WRITE(10,*) ' L2-norm on div of u =', norm_err(1)
383  WRITE(10,*) ' L2-norm of div of u err/norm =', norm_err(2)
384  WRITE(10,*) ' L2-norm of u =', norm_err(3)
385  WRITE(10,*) ' H1-norm of p =', norm_err(4)
386  WRITE(10,*) '########################################################'
387  END IF
388 
389  CASE(29)
390 
391  err = norm_sf(comm_one_d, 'div', vv_mesh, list_mode, un)
392  norm = norm_sf(comm_one_d, 'H1', vv_mesh, list_mode, un)
393  norm_err(1) = err/norm
394  norm = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, bn)
395  err = norm_sf(comm_one_d, 'div', h_mesh, list_mode, bn)
396  norm_err(2) = err/norm
397  norm_err(3) = norm
398  err= dot_product_sf(comm_one_d, h_mesh, list_mode, hn, bn)
399  norm_err(4) = 0.5*err
400  IF (rank==0) THEN
401  WRITE(10,*) '########################################################'
402  WRITE(10,*) ' L2-norm of div of u err/norm =', norm_err(1)
403  WRITE(10,*) ' L2-norm on div of B err/norm =', norm_err(2)
404  WRITE(10,*) ' L2-norm of of B =', norm_err(3)
405  WRITE(10,*) ' integral of 0.5*B.H =', norm_err(4)
406  WRITE(10,*) '########################################################'
407  END IF
408 
409  CASE(30)
410 
411  DO type=1,6
412  DO i=1,size(list_mode)
413  un_ex(:,TYPE,i) = vv_exact(TYPE,vv_mesh%rr,list_mode(i),time)
414  END DO
415  END DO
416  un_error = un - un_ex
417  norm_err(1) = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_error) / norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_ex)
418  DO type=1,2
419  DO i=1,size(list_mode)
420  pn_ex(:,TYPE,i) = pp_exact(TYPE,pp_mesh%rr,list_mode(i),time)
421  END DO
422  END DO
423  pn_error = pn - pn_ex
424  norm_err(2) = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, pn_error)
425  DO type=1,2
426  DO i=1,size(list_mode)
427  tempn_ex(:,TYPE,i) = temperature_exact(TYPE,temp_mesh%rr,list_mode(i),time)
428  END DO
429  END DO
430  tempn_error = tempn - tempn_ex
431  norm_err(3) = norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_error) / &
432  norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_ex) ! MODIFICATION: comm_one_d_temp instead of ns
433  norm_err(4) = norm_sf(comm_one_d_temp, 'H1', temp_mesh, list_mode, tempn_error) / &
434  norm_sf(comm_one_d_temp, 'H1', temp_mesh, list_mode, tempn_ex) ! MODIFICATION: comm_one_d_temp instead of ns
435  IF (rank==0) THEN
436  WRITE(10,*) '########################################################'
437  WRITE(10,*) ' L2-norm of error on u / L2-norm of u exact =', norm_err(1)
438  WRITE(10,*) ' L2-norm of error on p =', norm_err(2)
439  WRITE(10,*) ' L2-norm of error on T / L2-norm of T exact =', norm_err(3)
440  WRITE(10,*) ' H1-norm of error on T / H1-norm of T exact =', norm_err(4)
441  WRITE(10,*) '########################################################'
442  END IF
443 
444  CASE(31,32)
445 
446  DO type=1,6
447  DO i=1,size(list_mode)
448  un_ex(:,TYPE,i) = vv_exact(TYPE,vv_mesh%rr,list_mode(i),time)
449  END DO
450  END DO
451  un_error = un - un_ex
452  norm_err(1) = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_error) / norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_ex)
453  DO type=1,2
454  DO i=1,size(list_mode)
455  pn_ex(:,TYPE,i) = pp_exact(TYPE,pp_mesh%rr,list_mode(i),time)
456  END DO
457  END DO
458  pn_error = pn - pn_ex
459  norm_err(2) = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, pn_error) / norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, pn_ex)
460  DO type=1,2
461  DO i=1,size(list_mode)
462  tempn_ex(:,TYPE,i) = temperature_exact(TYPE,temp_mesh%rr,list_mode(i),time)
463  END DO
464  END DO
465  tempn_error = tempn - tempn_ex
466  norm_err(3) = norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_error) / &
467  norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_ex) ! MODIFICATION: comm_one_d_temp instead of ns
468  norm_err(4) = norm_sf(comm_one_d_temp, 'H1', temp_mesh, list_mode, tempn_error) / &
469  norm_sf(comm_one_d_temp, 'H1', temp_mesh, list_mode, tempn_ex) ! MODIFICATION: comm_one_d_temp instead of ns
470  IF (rank==0) THEN
471  WRITE(10,*) '########################################################'
472  WRITE(10,*) ' L2-norm of error on u / L2-norm of u exact =', norm_err(1)
473  WRITE(10,*) ' L2-norm of error on p / L2-norm of p exact =', norm_err(2)
474  WRITE(10,*) ' L2-norm of error on T / L2-norm of T exact =', norm_err(3)
475  WRITE(10,*) ' H1-norm of error on T / H1-norm of T exact =', norm_err(4)
476  WRITE(10,*) '########################################################'
477  END IF
478 
479  CASE(33)
480 
481  DO type=1,6
482  DO i=1,size(list_mode)
483  un_ex(:,TYPE,i) = vv_exact(TYPE,vv_mesh%rr,list_mode(i),time)
484  END DO
485  END DO
486  un_error = un - un_ex
487  norm_err(1) = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_error) / norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_ex)
488  DO type=1,2
489  DO i=1,size(list_mode)
490  pn_ex(:,TYPE,i) = pp_exact(TYPE,pp_mesh%rr,list_mode(i),time)
491  END DO
492  END DO
493  pn_error = pn - pn_ex
494  norm_err(2) = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, pn_error)
495  DO type=1,2
496  DO i=1,size(list_mode)
497  tempn_ex(:,TYPE,i) = temperature_exact(TYPE,temp_mesh%rr,list_mode(i),time)
498  END DO
499  END DO
500  tempn_error = tempn - tempn_ex
501  norm_err(3) = norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_error) / &
502  norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_ex) ! MODIFICATION: comm_one_d_temp instead of ns
503  DO type=1,6
504  DO i=1,size(list_mode)
505  hn_ex(:,TYPE,i) = hexact(h_mesh,TYPE,h_mesh%rr,list_mode(i),mu_h_field,time)
506  END DO
507  END DO
508  hn_error = hn - hn_ex
509  norm_err(4) = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn_error) / norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn_ex)
510  IF (rank==0) THEN
511  WRITE(10,*) '########################################################'
512  WRITE(10,*) ' L2-norm of error on u / L2-norm of u exact =', norm_err(1)
513  WRITE(10,*) ' L2-norm of error on p =', norm_err(2)
514  WRITE(10,*) ' L2-norm of error on T / L2-norm of T exact =', norm_err(3)
515  WRITE(10,*) ' L2-norm of error on H / L2-norm of H exact =', norm_err(4)
516  WRITE(10,*) '########################################################'
517  END IF
518 
519  CASE(34,37,38)
520 
521  DO type=1,6
522  DO i=1,size(list_mode)
523  un_ex(:,TYPE,i) = vv_exact(TYPE,vv_mesh%rr,list_mode(i),time)
524  END DO
525  END DO
526  un_error = un - un_ex
527  norm_err(1) = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_error) / norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_ex)
528  DO type=1,2
529  DO i=1,size(list_mode)
530  pn_ex(:,TYPE,i) = pp_exact(TYPE,pp_mesh%rr,list_mode(i),time)
531  END DO
532  END DO
533  pn_error = pn - pn_ex
534  norm_err(2) = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, pn_error) / norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, pn_ex)
535  DO type=1,2
536  DO i=1,size(list_mode)
537  tempn_ex(:,TYPE,i) = temperature_exact(TYPE,temp_mesh%rr,list_mode(i),time)
538  END DO
539  END DO
540  tempn_error = tempn - tempn_ex
541  norm_err(3) = norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_error) / &
542  norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_ex) ! MODIFICATION: comm_one_d_temp instead of ns
543  DO type=1,6
544  DO i=1,size(list_mode)
545  hn_ex(:,TYPE,i) = hexact(h_mesh,TYPE,h_mesh%rr,list_mode(i),mu_h_field,time)
546  END DO
547  END DO
548  hn_error = hn - hn_ex
549  norm_err(4) = norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn_error) / norm_sf(comm_one_d, 'L2', h_mesh, list_mode, hn_ex)
550  IF (rank==0) THEN
551  WRITE(10,*) '########################################################'
552  WRITE(10,*) ' L2-norm of error on u / L2-norm of u exact =', norm_err(1)
553  WRITE(10,*) ' L2-norm of error on p / L2-norm of p exact =', norm_err(2)
554  WRITE(10,*) ' L2-norm of error on T / L2-norm of T exact =', norm_err(3)
555  WRITE(10,*) ' L2-norm of error on H / L2-norm of H exact =', norm_err(4)
556  WRITE(10,*) '########################################################'
557  END IF
558 
559  CASE(35)
560 
561  DO type=1,6
562  DO i=1,size(list_mode)
563  un_ex(:,TYPE,i) = vv_exact(TYPE,vv_mesh%rr,list_mode(i),time)
564  END DO
565  END DO
566  un_error = un - un_ex
567  norm_err(1) = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_error)
568  DO type=1,2
569  DO i=1,size(list_mode)
570  pn_ex(:,TYPE,i) = pp_exact(TYPE,pp_mesh%rr,list_mode(i),time)
571  END DO
572  END DO
573  pn_error = pn - pn_ex
574  norm_err(2) = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, pn_error)
575  DO type=1,2
576  DO i=1,size(list_mode)
577  tempn_ex(:,TYPE,i) = temperature_exact(TYPE,temp_mesh%rr,list_mode(i),time)
578  END DO
579  END DO
580  tempn_error = tempn - tempn_ex
581  norm_err(3) = norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_error) / &
582  norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_ex)
583  norm_err(4) = norm_sf(comm_one_d_temp, 'H1', temp_mesh, list_mode, tempn_error) / &
584  norm_sf(comm_one_d_temp, 'H1', temp_mesh, list_mode, tempn_ex)
585  IF (rank==0) THEN
586  WRITE(10,*) '########################################################'
587  WRITE(10,*) ' L2-norm of error on u =', norm_err(1)
588  WRITE(10,*) ' L2-norm of error on p =', norm_err(2)
589  WRITE(10,*) ' L2-norm of error on T / L2-norm of T exact =', norm_err(3)
590  WRITE(10,*) ' H1-norm of error on T / H1-norm of T exact =', norm_err(4)
591  WRITE(10,*) '########################################################'
592  END IF
593 
594  CASE(36)
595 
596  DO type=1,6
597  DO i=1,size(list_mode)
598  un_ex(:,TYPE,i) = vv_exact(TYPE,vv_mesh%rr,list_mode(i),time)
599  END DO
600  END DO
601  un_error = un - un_ex
602  norm_err(1) = norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_error) / norm_sf(comm_one_d, 'L2', vv_mesh, list_mode, un_ex)
603  DO type=1,6
604  DO i=1,size(list_mode)
605  un_ex(:,TYPE,i) = vv_exact(TYPE,vv_mesh%rr,list_mode(i),time)
606  END DO
607  END DO
608  un_error = un - un_ex
609  norm_err(2) = norm_sf(comm_one_d, 'H1', vv_mesh, list_mode, un_error) / norm_sf(comm_one_d, 'H1', vv_mesh, list_mode, un_ex)
610  DO type=1,2
611  DO i=1,size(list_mode)
612  pn_ex(:,TYPE,i) = pp_exact(TYPE,pp_mesh%rr,list_mode(i),time)
613  END DO
614  END DO
615  pn_error = pn - pn_ex
616  norm_err(3) = norm_sf(comm_one_d, 'L2', pp_mesh, list_mode, pn_error)
617  DO type=1,2
618  DO i=1,size(list_mode)
619  tempn_ex(:,TYPE,i) = temperature_exact(TYPE,temp_mesh%rr,list_mode(i),time)
620  END DO
621  END DO
622  tempn_error = tempn - tempn_ex
623  norm_err(4) = norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_error) / &
624  norm_sf(comm_one_d_temp, 'L2', temp_mesh, list_mode, tempn_ex) ! MODIFICATION: comm_one_d_temp instead of ns
625  IF (rank==0) THEN
626  WRITE(10,*) '########################################################'
627  WRITE(10,*) ' L2-norm of error on u / L2-norm of u exact =', norm_err(1)
628  WRITE(10,*) ' H1-norm of error on u / H1-norm of u exact =', norm_err(2)
629  WRITE(10,*) ' L2-norm of error on p =', norm_err(3)
630  WRITE(10,*) ' L2-norm of error on T / L2-norm of T exact =', norm_err(4)
631  WRITE(10,*) '########################################################'
632  END IF
633 
634  CASE DEFAULT
635  CALL error_petsc(' BUG in post_proc_test: We should not be here')
636 
637  END SELECT
638 
639  IF (rank==0) THEN
640  ! Compare with references
641 !!$ IF (MAXVAL(ABS(inputs%norm_ref-norm_err)/inputs%norm_ref)<error_on_tests) THEN
642  IF (maxval(abs(inputs%norm_ref-norm_err))<1.d-8) THEN
643  WRITE(57,'(A,I2,A)') 'test #',inputs%numero_du_test_debug,' OK'
644  ELSE
645 !!$ WRITE(57,'(A,I2,2x,e15.7)') 'Problem with test #', inputs%numero_du_test_debug, &
646 !!$ MAXVAL(ABS(inputs%norm_ref-norm_err)/inputs%norm_ref)
647  WRITE(57,'(A,I2,2x,e15.7)') 'Problem with test #', inputs%numero_du_test_debug, &
648  maxval(abs(inputs%norm_ref-norm_err))
649  END IF
650  ! End compare
651  END IF
652 
653  END SUBROUTINE post_proc_test
654  !---------------------------------------------------------------------------
655 
656 END MODULE post_processing_debug
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
my_util
Definition: my_util.f90:1
post_processing_debug
Definition: post_processing_debug.f90:1
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:16
input_data::inputs
type(my_data), public inputs
Definition: read_my_data.f90:237
post_processing_debug::post_proc_test
subroutine, public post_proc_test(vv_mesh, pp_mesh, temp_mesh, H_mesh, phi_mesh, list_mode, un, pn, Hn, Bn, phin, tempn, level_setn, mu_H_field, time, m_max_c, comm_one_d, comm_one_d_ns, comm_one_d_temp)
Definition: post_processing_debug.f90:11
sft_parallele
Definition: fft_parallel.f90:4
tn_axi::angular_momentum
subroutine angular_momentum(mesh, list_mode, field, moments_out)
Definition: tn_axi.f90:257
tn_axi::norme_l2_champ_par
real(kind=8) function norme_l2_champ_par(communicator, mesh, list_mode, v)
Definition: tn_axi.f90:185
tn_axi
Definition: tn_axi.f90:5
tn_axi::norm_sf
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:40
subroutine_ns_with_u::moy
subroutine, public moy(communicator, mesh, p, RESLT)
Definition: sub_ns_with_velocity.f90:885
input_data
Definition: read_my_data.f90:233
tn_axi::dot_product_sf
real(kind=8) function dot_product_sf(communicator, mesh, list_mode, v, w)
Definition: tn_axi.f90:10
def_type_mesh
Definition: def_type_mesh.f90:4
tn_axi::norm_s_l1_zero_mode
real(kind=8) function norm_s_l1_zero_mode(communicator, mesh, list_mode, v)
Definition: tn_axi.f90:135
subroutine_ns_with_u
Definition: sub_ns_with_velocity.f90:5