From a0def9e9d81fd0c3806b220b02bfc5f50598c235 Mon Sep 17 00:00:00 2001
From: Joern Behrens <joern.behrens@uni-hamburg.de>
Date: Wed, 16 Feb 2022 18:28:02 +0100
Subject: [PATCH] added velocities to vtu output
---
.../ADV_semilagrange_stacked.F90 | 1 +
.../FLASH_parameters_stacked.f90 | 1 +
flash2d/src/options-sphere/IO_vtu_stacked.F90 | 37 ++++++++++++++++++-
.../src/options-sphere/IO_vtuplot_stacked.F90 | 32 ++++++++++++++--
.../src/options-sphere/SLM_simple_stacked.f90 | 14 ++++---
5 files changed, 74 insertions(+), 11 deletions(-)
diff --git a/flash2d/src/options-sphere/ADV_semilagrange_stacked.F90 b/flash2d/src/options-sphere/ADV_semilagrange_stacked.F90
index 2972464..e804907 100644
--- a/flash2d/src/options-sphere/ADV_semilagrange_stacked.F90
+++ b/flash2d/src/options-sphere/ADV_semilagrange_stacked.F90
@@ -523,6 +523,7 @@
i_stcnumlayers = p_param%tst%tst_int(1,1)
i_stctracer = grid_registerfemvar(1,i_length=i_stcnumlayers)
+ i_stcuv = grid_registerfemvar(1,i_length=(i_stcnumlayers*3))
allocate(r_stclayerpressure(i_stcnumlayers), stat=i_alct)
not_allc: IF(i_alct /= 0) THEN
CALL grid_error(c_error='[slm_initialize]: Could not allocate layers pressure array')
diff --git a/flash2d/src/options-sphere/FLASH_parameters_stacked.f90 b/flash2d/src/options-sphere/FLASH_parameters_stacked.f90
index 9dd8c8c..d58e3c2 100644
--- a/flash2d/src/options-sphere/FLASH_parameters_stacked.f90
+++ b/flash2d/src/options-sphere/FLASH_parameters_stacked.f90
@@ -113,6 +113,7 @@
!---------- convenience pointers to stacked layers
INTEGER (KIND = GRID_SI) :: i_stctracer
+ INTEGER (KIND = GRID_SI) :: i_stcuv
INTEGER (KIND = GRID_SI) :: i_stcnumlayers
!---------- array for storing stacked layer pressure levels
diff --git a/flash2d/src/options-sphere/IO_vtu_stacked.F90 b/flash2d/src/options-sphere/IO_vtu_stacked.F90
index d589285..08e0df5 100644
--- a/flash2d/src/options-sphere/IO_vtu_stacked.F90
+++ b/flash2d/src/options-sphere/IO_vtu_stacked.F90
@@ -135,6 +135,9 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
REAL ( KIND = GRID_SR ) :: r_heightfactor, r_aux
TYPE(t_vtu_data), DIMENSION(6) :: p_stcnodedata
REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_stcval
+ REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_stcuwnd
+ REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_stcvwnd
+ REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_stcwwnd
INTEGER( KIND = GRID_SI) :: i_stcbeg
INTEGER( KIND = GRID_SI) :: i_stcend
CHARACTER (len=12) :: c_layer
@@ -406,7 +409,10 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
CALL write_vtu_data(i_fhandle, p_nodedata(i_cnt))
END DO
#elif defined VTU_OUTPUTSTC
- ALLOCATE(r_stcval(1,i_numberofvtupoints), stat=i_alct)
+ ALLOCATE(r_stcval(1,i_numberofvtupoints), &
+ r_stcuwnd(1,i_numberofvtupoints), &
+ r_stcvwnd(1,i_numberofvtupoints), &
+ r_stcwwnd(1,i_numberofvtupoints), stat=i_alct)
IF(i_alct /= 0) THEN
CALL grid_error(c_error='[plot_vtu]: could not allocate aux. stc values arrays')
END IF
@@ -417,13 +423,40 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
i_stcbeg = p_mesh%i_nnumber*(i_hcnt-1)+1
i_stcend = p_mesh%i_nnumber*i_hcnt
r_stcval(1:1,i_stcbeg:i_stcend) = p_nodedata(i_cnt)%p_vdata(1:1,:)
+ ELSE IF (p_nodedata(i_cnt)%c_name(1:6) == ' u-wnd') THEN
+ c_layer = p_nodedata(i_cnt)%c_name(7:22)
+ READ(c_layer,*) i_hcnt
+ i_stcbeg = p_mesh%i_nnumber*(i_hcnt-1)+1
+ i_stcend = p_mesh%i_nnumber*i_hcnt
+ r_stcuwnd(1:1,i_stcbeg:i_stcend) = p_nodedata(i_cnt)%p_vdata(1:1,:)
+ ELSE IF (p_nodedata(i_cnt)%c_name(1:6) == ' v-wnd') THEN
+ c_layer = p_nodedata(i_cnt)%c_name(7:22)
+ READ(c_layer,*) i_hcnt
+ i_stcbeg = p_mesh%i_nnumber*(i_hcnt-1)+1
+ i_stcend = p_mesh%i_nnumber*i_hcnt
+ r_stcvwnd(1:1,i_stcbeg:i_stcend) = p_nodedata(i_cnt)%p_vdata(1:1,:)
+ ELSE IF (p_nodedata(i_cnt)%c_name(1:6) == ' w-wnd') THEN
+ c_layer = p_nodedata(i_cnt)%c_name(7:22)
+ READ(c_layer,*) i_hcnt
+ i_stcbeg = p_mesh%i_nnumber*(i_hcnt-1)+1
+ i_stcend = p_mesh%i_nnumber*i_hcnt
+ r_stcwwnd(1:1,i_stcbeg:i_stcend) = p_nodedata(i_cnt)%p_vdata(1:1,:)
END IF
END DO
p_stcnodedata(1)%c_name = 'stc-tracer'
p_stcnodedata(1)%i_size = 1
p_stcnodedata(1)%p_vdata=>r_stcval
CALL write_vtu_data(i_fhandle, p_stcnodedata(1))
- DEALLOCATE(r_stcval)
+ p_stcnodedata(1)%c_name = 'stc-u-wind'
+ p_stcnodedata(1)%p_vdata=>r_stcuwnd
+ CALL write_vtu_data(i_fhandle, p_stcnodedata(1))
+ p_stcnodedata(1)%c_name = 'stc-v-wind'
+ p_stcnodedata(1)%p_vdata=>r_stcvwnd
+ CALL write_vtu_data(i_fhandle, p_stcnodedata(1))
+ p_stcnodedata(1)%c_name = 'stc-w-wind'
+ p_stcnodedata(1)%p_vdata=>r_stcwwnd
+ CALL write_vtu_data(i_fhandle, p_stcnodedata(1))
+ DEALLOCATE(r_stcval, r_stcuwnd, r_stcvwnd, r_stcwwnd)
#endif
END IF
diff --git a/flash2d/src/options-sphere/IO_vtuplot_stacked.F90 b/flash2d/src/options-sphere/IO_vtuplot_stacked.F90
index c1ea6cc..52dec7a 100644
--- a/flash2d/src/options-sphere/IO_vtuplot_stacked.F90
+++ b/flash2d/src/options-sphere/IO_vtuplot_stacked.F90
@@ -55,10 +55,11 @@ MODULE IO_vtuplot
INTEGER, PARAMETER :: i_vallen=5
INTEGER (KIND = GRID_SI), DIMENSION(:), ALLOCATABLE :: i_valind
INTEGER (KIND = GRID_SI) :: i_trcvals, i_lay, i_nodnum
+ INTEGER (KIND = GRID_SI) :: i_velovals, i_vcnt
TYPE(t_vtu_data), DIMENSION(2) :: celldata
- ! CAUTION: This is a fixed dimension, but will only work for max 16 layers...
- TYPE(t_vtu_data), DIMENSION(20) :: nodedata
+ ! CAUTION: This is a fixed dimension, but will only work for max 11 layers...
+ TYPE(t_vtu_data), DIMENSION(50) :: nodedata
!---------- file handling (open)
@@ -101,13 +102,19 @@ MODULE IO_vtuplot
i_valind= (/GRID_ucomp, GRID_vcomp, GRID_phi, GRID_zeta, GRID_tracer, &
((i_stctracer-1+i_lay), i_lay=1,i_stcnumlayers) /)
CALL grid_getinfo(p_handle, r_nodevalues= r_val, i_arraypoint=i_valind)
+
+!---------- determine number of velocity values from layers and allocate accordingly
+
+ i_velovals = GRID_dimension * i_stcnumlayers
!---------- assign nodal values, but before allocate aux. arrays
- ALLOCATE(r_velo(2,i_nnum), stat=i_alct)
+ ALLOCATE(r_velo(i_velovals,i_nnum), stat=i_alct)
IF(i_alct /= 0) THEN
CALL grid_error(c_error='[generate_vtu]: could not allocate aux. velo arrays')
END IF
+ i_valind= (/ ((i_stcuv-1+i_lay), i_lay=1,i_velovals) /)
+ CALL grid_getinfo(p_handle, r_nodevalues= r_velo, i_arraypoint=i_valind)
!---------- generate 2D velocity field
@@ -141,6 +148,25 @@ MODULE IO_vtuplot
END DO
i_nodnum= i_stcnumlayers+ 5
+!---------- encode velocity fields for layers
+
+ DO i_lay=1,i_stcnumlayers
+ i_vcnt= (i_lay-1)*GRID_dimension+1+i_nodnum
+ write(c_tmp,*) 'u-wnd',i_lay
+ nodedata(i_lay+5)%c_name = c_tmp
+ nodedata(i_lay+5)%i_size = 1
+ nodedata(i_lay+5)%p_vdata => r_velo(i_vcnt:i_vcnt,:)
+ write(c_tmp,*) 'v-wnd',i_lay
+ nodedata(i_lay+5)%c_name = c_tmp
+ nodedata(i_lay+5)%i_size = 1
+ nodedata(i_lay+5)%p_vdata => r_velo(i_vcnt+1:i_vcnt+1,:)
+ write(c_tmp,*) 'w-wnd',i_lay
+ nodedata(i_lay+5)%c_name = c_tmp
+ nodedata(i_lay+5)%i_size = 1
+ nodedata(i_lay+5)%p_vdata => r_velo(i_vcnt+2:i_vcnt+2,:)
+ END DO
+ i_nodnum= i_nodnum+ GRID_dimension*i_stcnumlayers
+
!---------- extract cells grid data
i_tnum= p_handle%i_enumfine
diff --git a/flash2d/src/options-sphere/SLM_simple_stacked.f90 b/flash2d/src/options-sphere/SLM_simple_stacked.f90
index 4c42f69..cfad90f 100644
--- a/flash2d/src/options-sphere/SLM_simple_stacked.f90
+++ b/flash2d/src/options-sphere/SLM_simple_stacked.f90
@@ -145,6 +145,8 @@
REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE :: r_newvl
REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_alpha
INTEGER (KIND = GRID_SI) :: i_alct, i_layerid, i_nd
+ INTEGER (KIND = GRID_SI), DIMENSION(GRID_dimension) :: i_uvpointer
+ INTEGER (KIND = GRID_SI) :: i_start, i_cnt
!---------- evaluate layer information, if given
@@ -195,15 +197,15 @@
!-SLM--------- put alpha values to u and v field entries
- IF(i_layerid == 1) THEN
+ i_start= i_layerid-1
+ i_uvpointer = (/ (i_start+i_cnt, i_cnt=1,GRID_dimension) /)
r_alpha= -r_alpha
IF(present(i_newsdepth)) THEN
- CALL grid_putinfo(p_ghand(i_timeplus), r_nodevalues=r_alpha(1:2,:), &
- i_newsdepth=i_newsdepth, i_arraypoint=(/ GRID_ucomp, GRID_vcomp /))
+ CALL grid_putinfo(p_ghand(i_timeplus), r_nodevalues=r_alpha(1:GRID_dimension,:), &
+ i_newsdepth=i_newsdepth, i_arraypoint=i_uvpointer)
ELSE
- CALL grid_putinfo(p_ghand(i_timeplus), r_nodevalues=r_alpha(1:2,:), &
- i_arraypoint=(/ GRID_ucomp, GRID_vcomp /))
- END IF
+ CALL grid_putinfo(p_ghand(i_timeplus), r_nodevalues=r_alpha(1:GRID_dimension,:), &
+ i_arraypoint=i_uvpointer)
END IF
!-SLM--------- deallocate work arrays
--
GitLab