Write low-Mach compatibility and source-time-level diagnostics.
The current divergence_source may be advanced by the energy/thermo density
sync after the projection has already used the previous source level. The
projection_divergence_source field stores the source that was actually
consumed by the latest projection RHS and outlet balancing. Comparing both
residuals distinguishes pressure/projection error from source time-level
mismatch.
CSV columns with _current compare against fields%divergence_source
after the energy/thermo update; columns with _projection compare
against the source copied before projection.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(mesh_t), | intent(in) | :: | mesh | |||
| type(flow_mpi_t), | intent(in) | :: | flow | |||
| type(case_params_t), | intent(in) | :: | params | |||
| type(flow_fields_t), | intent(in) | :: | fields | |||
| integer, | intent(in) | :: | step | |||
| real(kind=rk), | intent(in) | :: | time |
subroutine write_variable_density_compatibility_diagnostics(mesh, flow, params, fields, step, time) type(mesh_t), intent(in) :: mesh type(flow_mpi_t), intent(in) :: flow type(case_params_t), intent(in) :: params type(flow_fields_t), intent(in) :: fields integer, intent(in) :: step real(rk), intent(in) :: time integer :: c, f, lf, fid, owner, nb, ierr, unit_id logical, save :: compatibility_diag_initialized = .false. character(len=1024) :: filename real(rk) :: vol, s_c, s_proj_c, div_c, res_current, res_projection, s_delta, flux_out real(rk) :: local_integral_S_dV, integral_S_dV real(rk) :: local_integral_S_projection_dV, integral_S_projection_dV real(rk) :: local_net_boundary_volume_flux, net_boundary_volume_flux real(rk) :: compatibility_error, compatibility_error_projection real(rk) :: local_volume, volume_total, mean_S, mean_S_projection real(rk) :: local_S_l2_num, S_l2_num, S_l2 real(rk) :: local_S_projection_l2_num, S_projection_l2_num, S_projection_l2 real(rk) :: local_S_abs_max, S_abs_max real(rk) :: local_S_projection_abs_max, S_projection_abs_max real(rk) :: local_res_max, res_max real(rk) :: local_res_projection_max, res_projection_max real(rk) :: local_res_l2_num, res_l2_num, res_l2 real(rk) :: local_res_projection_l2_num, res_projection_l2_num, res_projection_l2 real(rk) :: local_S_delta_max, S_delta_max real(rk) :: local_S_delta_l2_num, S_delta_l2_num, S_delta_l2 real(rk) :: rel_res_max, rel_res_l2 real(rk) :: rel_res_projection_max, rel_res_projection_l2 real(rk) :: denom real(rk) :: mpi_reduce_send(1), mpi_reduce_recv(1) if (.not. params%enable_variable_density) return if (.not. params%write_diagnostics) return local_integral_S_dV = 0.0_rk local_integral_S_projection_dV = 0.0_rk local_net_boundary_volume_flux = 0.0_rk local_volume = 0.0_rk local_S_l2_num = 0.0_rk local_S_projection_l2_num = 0.0_rk local_S_abs_max = 0.0_rk local_S_projection_abs_max = 0.0_rk local_res_max = 0.0_rk local_res_projection_max = 0.0_rk local_res_l2_num = 0.0_rk local_res_projection_l2_num = 0.0_rk local_S_delta_max = 0.0_rk local_S_delta_l2_num = 0.0_rk do c = flow%first_cell, flow%last_cell vol = mesh%cells(c)%volume s_c = 0.0_rk if (allocated(fields%divergence_source)) s_c = fields%divergence_source(c) s_proj_c = s_c if (allocated(fields%projection_divergence_source)) s_proj_c = fields%projection_divergence_source(c) div_c = 0.0_rk if (allocated(fields%face_flux)) then do lf = 1, mesh%ncell_faces(c) fid = mesh%cell_faces(lf, c) if (mesh%faces(fid)%owner == c) then flux_out = fields%face_flux(fid) else flux_out = -fields%face_flux(fid) end if div_c = div_c + flux_out / vol end do end if res_current = div_c - s_c res_projection = div_c - s_proj_c s_delta = s_c - s_proj_c local_integral_S_dV = local_integral_S_dV + s_c * vol local_integral_S_projection_dV = local_integral_S_projection_dV + s_proj_c * vol local_volume = local_volume + vol local_S_l2_num = local_S_l2_num + s_c * s_c * vol local_S_projection_l2_num = local_S_projection_l2_num + s_proj_c * s_proj_c * vol local_S_abs_max = max(local_S_abs_max, abs(s_c)) local_S_projection_abs_max = max(local_S_projection_abs_max, abs(s_proj_c)) local_res_max = max(local_res_max, abs(res_current)) local_res_projection_max = max(local_res_projection_max, abs(res_projection)) local_res_l2_num = local_res_l2_num + res_current * res_current * vol local_res_projection_l2_num = local_res_projection_l2_num + res_projection * res_projection * vol local_S_delta_max = max(local_S_delta_max, abs(s_delta)) local_S_delta_l2_num = local_S_delta_l2_num + s_delta * s_delta * vol end do if (allocated(fields%face_flux)) then do f = 1, mesh%nfaces owner = mesh%faces(f)%owner nb = mesh%faces(f)%neighbor if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor if (nb <= 0) then if (owner >= flow%first_cell .and. owner <= flow%last_cell) then local_net_boundary_volume_flux = local_net_boundary_volume_flux + fields%face_flux(f) end if end if end do end if mpi_reduce_send(1) = local_integral_S_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_S_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility integral S dV') mpi_reduce_send(1) = local_integral_S_projection_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_S_projection_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility projection integral S dV') mpi_reduce_send(1) = local_net_boundary_volume_flux call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) net_boundary_volume_flux = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility boundary flux') mpi_reduce_send(1) = local_volume call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) volume_total = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility volume') mpi_reduce_send(1) = local_S_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) S_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility S l2') mpi_reduce_send(1) = local_S_projection_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) S_projection_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility projection S l2') mpi_reduce_send(1) = local_S_abs_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) S_abs_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility S max') mpi_reduce_send(1) = local_S_projection_abs_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) S_projection_abs_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility projection S max') mpi_reduce_send(1) = local_res_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) res_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility residual max') mpi_reduce_send(1) = local_res_projection_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) res_projection_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility projection residual max') mpi_reduce_send(1) = local_res_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) res_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility residual l2') mpi_reduce_send(1) = local_res_projection_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) res_projection_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility projection residual l2') mpi_reduce_send(1) = local_S_delta_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) S_delta_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility source delta max') mpi_reduce_send(1) = local_S_delta_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) S_delta_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing compatibility source delta l2') compatibility_error = net_boundary_volume_flux - integral_S_dV compatibility_error_projection = net_boundary_volume_flux - integral_S_projection_dV if (volume_total > 0.0_rk) then mean_S = integral_S_dV / volume_total mean_S_projection = integral_S_projection_dV / volume_total S_l2 = sqrt(S_l2_num / volume_total) S_projection_l2 = sqrt(S_projection_l2_num / volume_total) res_l2 = sqrt(res_l2_num / volume_total) res_projection_l2 = sqrt(res_projection_l2_num / volume_total) S_delta_l2 = sqrt(S_delta_l2_num / volume_total) else mean_S = 0.0_rk mean_S_projection = 0.0_rk S_l2 = 0.0_rk S_projection_l2 = 0.0_rk res_l2 = 0.0_rk res_projection_l2 = 0.0_rk S_delta_l2 = 0.0_rk end if denom = max(S_abs_max, tiny(1.0_rk)) rel_res_max = res_max / denom denom = max(S_l2, tiny(1.0_rk)) rel_res_l2 = res_l2 / denom denom = max(S_projection_abs_max, tiny(1.0_rk)) rel_res_projection_max = res_projection_max / denom denom = max(S_projection_l2, tiny(1.0_rk)) rel_res_projection_l2 = res_projection_l2 / denom if (flow%rank == 0) then filename = trim(params%output_dir) // '/diagnostics/variable_density_compatibility.csv' if (.not. compatibility_diag_initialized) then open(newunit=unit_id, file=trim(filename), status='replace', action='write') compatibility_diag_initialized = .true. write(unit_id,'(a)') 'step,time,integral_S_current_dV,net_boundary_volume_flux,' // & 'net_boundary_volume_flux_minus_integral_S_current_dV,volume_total,mean_S_current,' // & 'S_current_l2,S_current_abs_max,divu_minus_S_current_max,divu_minus_S_current_l2,' // & 'relative_divu_minus_S_current_max,relative_divu_minus_S_current_l2,' // & 'integral_S_projection_dV,' // & 'net_boundary_volume_flux_minus_integral_S_projection_dV,' // & 'mean_S_projection,S_projection_l2,S_projection_abs_max,' // & 'divu_minus_S_projection_max,divu_minus_S_projection_l2,' // & 'relative_divu_minus_S_projection_max,' // & 'relative_divu_minus_S_projection_l2,' // & 'S_current_minus_S_projection_max,S_current_minus_S_projection_l2' else open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write') end if write(unit_id,'(i0,23(",",ES26.16E4))') step, time, integral_S_dV, net_boundary_volume_flux, & compatibility_error, volume_total, mean_S, S_l2, S_abs_max, & res_max, res_l2, rel_res_max, rel_res_l2, & integral_S_projection_dV, compatibility_error_projection, & mean_S_projection, S_projection_l2, S_projection_abs_max, & res_projection_max, res_projection_l2, & rel_res_projection_max, rel_res_projection_l2, & S_delta_max, S_delta_l2 close(unit_id) end if end subroutine write_variable_density_compatibility_diagnostics