diagnostics/m_diagnostics_binary.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Module for binary output | ||
| 2 | module m_diagnostics_binary | ||
| 3 | use m_common, only: wp | ||
| 4 | implicit none | ||
| 5 | |||
| 6 | private | ||
| 7 | public :: Binary2D | ||
| 8 | |||
| 9 | !> @brief Type for binary output of 2D m-forms | ||
| 10 | type :: Binary2D | ||
| 11 | integer :: current_file_index = 0 | ||
| 12 | character(len=256) :: filename_base | ||
| 13 | logical :: is_initialized = .false. | ||
| 14 | |||
| 15 | real(wp) :: xlims(2), ylims(2) | ||
| 16 | contains | ||
| 17 | procedure :: init => init_binary | ||
| 18 | procedure :: write => write_binary | ||
| 19 | end type Binary2D | ||
| 20 | contains | ||
| 21 | !> @brief Initialize the binary output | ||
| 22 | !> | ||
| 23 | !> @param[inout] this The Binary2D object | ||
| 24 | !> @param[in] filename Base filename for output files | ||
| 25 | !> @param[in] xlims Array of size 2 specifying the x-limits in physical space | ||
| 26 | !> @param[in] ylims Array of size 2 specifying the y-limits in physical space | ||
| 27 | ✗ | subroutine init_binary(this, filename, xlims, ylims) | |
| 28 | implicit none | ||
| 29 | |||
| 30 | class(Binary2D), intent(inout) :: this | ||
| 31 | character(len=*), intent(in) :: filename | ||
| 32 | real(wp), intent(in) :: xlims(2), ylims(2) | ||
| 33 | |||
| 34 | ✗ | this%filename_base = filename | |
| 35 | ✗ | this%current_file_index = 0 | |
| 36 | ✗ | this%is_initialized = .true. | |
| 37 | |||
| 38 | ✗ | this%xlims = xlims | |
| 39 | ✗ | this%ylims = ylims | |
| 40 | ✗ | end subroutine init_binary | |
| 41 | |||
| 42 | !> @brief Write the m-form values to a binary file | ||
| 43 | !> | ||
| 44 | !> @param[inout] this The Binary2D object | ||
| 45 | !> @param[in] mfun The m-form function to evaluate and write | ||
| 46 | !> @param[in] coord_transform _(optional)_ The coordinate transformation associated with the m-form | ||
| 47 | !> @param[in] refine _(optional)_ integer specifying refinement factor for output grid | ||
| 48 | ✗ | subroutine write_binary(this, mfun, coord_transform, refine) | |
| 49 | use, intrinsic :: ieee_arithmetic, only: ieee_quiet_nan, ieee_value | ||
| 50 | use m_mform_basis, only: MFormFun, evaluate | ||
| 51 | use m_coord_transform_abstract, only: CoordTransformAbstract | ||
| 52 | use m_coord_transform_interface, only: inverse_transform | ||
| 53 | implicit none | ||
| 54 | |||
| 55 | class(Binary2D), intent(inout) :: this | ||
| 56 | type(MFormFun), intent(in) :: mfun | ||
| 57 | class(CoordTransformAbstract), optional, intent(in) :: coord_transform | ||
| 58 | integer, optional, intent(in) :: refine | ||
| 59 | |||
| 60 | character(len=256) :: fname | ||
| 61 | integer :: unit, ierr, refine_, nx, ny, nr_blocks, blk, i, j | ||
| 62 | real(wp) :: x0, x1, y0, y1, dx, dy, x, y, val(3), xp, yp, zp | ||
| 63 | |||
| 64 | ✗ | if (.not. this%is_initialized) error stop "Binary2D::write: object not initialized." | |
| 65 | |||
| 66 | ✗ | write (fname, '(A,I5.5,A)') trim(this%filename_base), this%current_file_index, '.dat' | |
| 67 | |||
| 68 | refine_ = 0 | ||
| 69 | ✗ | if (present(refine)) refine_ = refine | |
| 70 | |||
| 71 | ! Get number of intervals from the space | ||
| 72 | ✗ | nx = mfun%space%tp_spaces(1)%spaces(1)%nr_intervals * (1 + refine_) + 1 | |
| 73 | ✗ | ny = mfun%space%tp_spaces(1)%spaces(2)%nr_intervals * (1 + refine_) + 1 | |
| 74 | |||
| 75 | ✗ | x0 = this%xlims(1) | |
| 76 | ✗ | x1 = this%xlims(2) | |
| 77 | ✗ | y0 = this%ylims(1) | |
| 78 | ✗ | y1 = this%ylims(2) | |
| 79 | |||
| 80 | ✗ | dx = (x1 - x0) / real(nx - 1, wp) | |
| 81 | ✗ | dy = (y1 - y0) / real(ny - 1, wp) | |
| 82 | |||
| 83 | ✗ | open (newunit=unit, file=fname, form='unformatted', access='stream', status='replace', action='write', iostat=ierr) | |
| 84 | ✗ | if (ierr /= 0) error stop 'Could not open file for binary grid output.' | |
| 85 | |||
| 86 | ✗ | write (unit) mfun%space%m | |
| 87 | |||
| 88 | ✗ | if (mfun%space%m == 1) then | |
| 89 | nr_blocks = 2 | ||
| 90 | else ! m = 0, 2 | ||
| 91 | nr_blocks = 1 | ||
| 92 | end if | ||
| 93 | |||
| 94 | ! Write header: number of points in each direction | ||
| 95 | ✗ | write (unit) nx, ny | |
| 96 | |||
| 97 | ✗ | write (unit) x0, x1, y0, y1 | |
| 98 | |||
| 99 | ! Write m-form values at each grid point (row-major order) | ||
| 100 | ✗ | do j = 0, ny - 1 | |
| 101 | ✗ | y = y0 + j * dy | |
| 102 | ✗ | do i = 0, nx - 1 | |
| 103 | ✗ | x = x0 + i * dx | |
| 104 | |||
| 105 | ! Shift to avoid pushforward problems at the axis | ||
| 106 | ✗ | if (present(coord_transform)) then | |
| 107 | ✗ | xp = inverse_transform(coord_transform, 1, x, y, 0.0_wp) | |
| 108 | ✗ | if (coord_transform%has_polar_xy_singularity .and. abs(xp) <= epsilon(100._wp)) xp = xp + sqrt(epsilon(1.0_wp)) | |
| 109 | ✗ | yp = inverse_transform(coord_transform, 2, x, y, 0.0_wp) | |
| 110 | ✗ | zp = inverse_transform(coord_transform, 3, x, y, 0.0_wp) | |
| 111 | else | ||
| 112 | ✗ | xp = x | |
| 113 | ✗ | yp = y | |
| 114 | ✗ | zp = 0.0_wp | |
| 115 | end if | ||
| 116 | |||
| 117 | ✗ | if (mfun%space%tp_spaces(1)%spaces(1)%is_periodic) then | |
| 118 | ✗ | xp = modulo(xp, 1.0_wp) | |
| 119 | end if | ||
| 120 | ✗ | if (mfun%space%tp_spaces(1)%spaces(2)%is_periodic) then | |
| 121 | ✗ | yp = modulo(yp, 1.0_wp) | |
| 122 | end if | ||
| 123 | |||
| 124 | ✗ | if (xp >= 0.0_wp .and. xp <= 1.0_wp .and. yp >= 0.0_wp .and. yp <= 1.0_wp) then | |
| 125 | ✗ | if (mfun%space%m == 0 .or. mfun%space%m == 3) then | |
| 126 | ✗ | val(1) = evaluate(mfun, xp, yp, zp, coord_transform=coord_transform) | |
| 127 | ✗ | elseif (mfun%space%m == 1) then | |
| 128 | ✗ | do blk = 1, nr_blocks | |
| 129 | ✗ | val(blk) = evaluate(mfun, blk, xp, yp, zp, coord_transform=coord_transform) | |
| 130 | end do | ||
| 131 | else ! m = 2 | ||
| 132 | ✗ | val(1) = evaluate(mfun, 3, xp, yp, zp, coord_transform=coord_transform) | |
| 133 | end if | ||
| 134 | else | ||
| 135 | ✗ | val = ieee_value(val, ieee_quiet_nan) | |
| 136 | end if | ||
| 137 | |||
| 138 | ✗ | write (unit) val(1:nr_blocks) | |
| 139 | end do | ||
| 140 | end do | ||
| 141 | |||
| 142 | ✗ | close (unit) | |
| 143 | |||
| 144 | ✗ | this%current_file_index = this%current_file_index + 1 | |
| 145 | |||
| 146 | ✗ | end subroutine write_binary | |
| 147 | ✗ | end module m_diagnostics_binary | |
| 148 |