Write variable-density transport/conservation diagnostics.
Boundary mass flux is positive outward, so global conservative closure is:
d/dt integral(rho dV) + net_boundary_mass_flux = 0
The time derivative here is output-to-output, not per-step. This keeps the diagnostic independent of the integrator state. Write variable-density mass/transport conservation diagnostics with explicit density and mass-flux time levels.
Sign convention: positive boundary flux is outward from the owner cell/domain.
For a conservative domain mass balance: dM/dt + net_boundary_mass_flux = 0
This routine reports several versions of the boundary mass flux so the diagnostics can distinguish a true conservative transport error from a current-density/projection-density time-level mismatch.
| 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 | |||
| type(energy_fields_t), | intent(in) | :: | energy | |||
| type(transport_properties_t), | intent(in) | :: | transport | |||
| integer, | intent(in) | :: | step | |||
| real(kind=rk), | intent(in) | :: | time |
subroutine write_variable_density_transport_conservation_diagnostics(mesh, flow, params, fields, energy, transport, 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 type(energy_fields_t), intent(in) :: energy type(transport_properties_t), intent(in) :: transport integer, intent(in) :: step real(rk), intent(in) :: time integer :: c, f, owner, nb, ierr, unit_id logical, save :: transport_diag_initialized = .false. character(len=1024) :: filename real(rk), save :: prev_time = 0.0_rk real(rk), save :: prev_mass_current = 0.0_rk real(rk), save :: prev_mass_projection = 0.0_rk real(rk), save :: prev_rho_h_current = 0.0_rk real(rk) :: vol, rho_c, rho_p, h_c, T_c real(rk) :: face_flux, mass_flux_stored real(rk) :: dt_prev real(rk) :: delta_mass_current, mass_rate_current real(rk) :: delta_mass_projection, mass_rate_projection real(rk) :: delta_rho_h_current, rho_h_rate_current real(rk) :: defect_current_stored, defect_current_current_rho real(rk) :: defect_projection_projection_rho real(rk) :: rel_defect_current_stored, rel_defect_current_current_rho real(rk) :: rel_defect_projection_projection_rho real(rk) :: denom real(rk) :: sum_send(12), sum_recv(12) real(rk) :: min_send(4), min_recv(4) real(rk) :: max_send(5), max_recv(5) real(rk) :: local_volume real(rk) :: local_mass_current, local_mass_projection real(rk) :: local_net_boundary_volume_flux real(rk) :: local_net_boundary_mass_flux_stored real(rk) :: local_net_boundary_mass_flux_current_rho real(rk) :: local_net_boundary_mass_flux_projection_rho real(rk) :: local_rho_current_l2_num, local_rho_projection_l2_num real(rk) :: local_rho_diff_l2_num real(rk) :: local_rho_h_current, local_rho_h_projection real(rk) :: volume_total real(rk) :: mass_current, mass_projection real(rk) :: current_minus_projection_mass real(rk) :: net_boundary_volume_flux real(rk) :: net_boundary_mass_flux_stored real(rk) :: net_boundary_mass_flux_current_rho real(rk) :: net_boundary_mass_flux_projection_rho real(rk) :: stored_minus_current_rho_boundary_mass_flux real(rk) :: stored_minus_projection_rho_boundary_mass_flux real(rk) :: rho_current_min, rho_current_max, rho_current_mean, rho_current_l2 real(rk) :: rho_projection_min, rho_projection_max, rho_projection_mean, rho_projection_l2 real(rk) :: rho_diff_max, rho_diff_l2 real(rk) :: h_min, h_max, T_min, T_max real(rk) :: rho_h_current, rho_h_projection if (.not. params%enable_variable_density) return if (.not. params%write_diagnostics) return if (.not. allocated(transport%rho)) return local_volume = 0.0_rk local_mass_current = 0.0_rk local_mass_projection = 0.0_rk local_net_boundary_volume_flux = 0.0_rk local_net_boundary_mass_flux_stored = 0.0_rk local_net_boundary_mass_flux_current_rho = 0.0_rk local_net_boundary_mass_flux_projection_rho = 0.0_rk local_rho_current_l2_num = 0.0_rk local_rho_projection_l2_num = 0.0_rk local_rho_diff_l2_num = 0.0_rk local_rho_h_current = 0.0_rk local_rho_h_projection = 0.0_rk rho_current_min = huge(1.0_rk) rho_projection_min = huge(1.0_rk) h_min = huge(1.0_rk) T_min = huge(1.0_rk) rho_current_max = -huge(1.0_rk) rho_projection_max = -huge(1.0_rk) rho_diff_max = 0.0_rk h_max = -huge(1.0_rk) T_max = -huge(1.0_rk) do c = flow%first_cell, flow%last_cell vol = mesh%cells(c)%volume rho_c = transport%rho(c) rho_p = rho_c if (allocated(fields%projection_rho)) then if (fields%projection_rho(c) > 0.0_rk) rho_p = fields%projection_rho(c) end if h_c = 0.0_rk if (allocated(energy%h)) h_c = energy%h(c) T_c = 0.0_rk if (allocated(energy%T)) T_c = energy%T(c) local_volume = local_volume + vol local_mass_current = local_mass_current + rho_c * vol local_mass_projection = local_mass_projection + rho_p * vol local_rho_current_l2_num = local_rho_current_l2_num + rho_c * rho_c * vol local_rho_projection_l2_num = local_rho_projection_l2_num + rho_p * rho_p * vol local_rho_diff_l2_num = local_rho_diff_l2_num + (rho_c - rho_p) * (rho_c - rho_p) * vol local_rho_h_current = local_rho_h_current + rho_c * h_c * vol local_rho_h_projection = local_rho_h_projection + rho_p * h_c * vol rho_current_min = min(rho_current_min, rho_c) rho_current_max = max(rho_current_max, rho_c) rho_projection_min = min(rho_projection_min, rho_p) rho_projection_max = max(rho_projection_max, rho_p) rho_diff_max = max(rho_diff_max, abs(rho_c - rho_p)) h_min = min(h_min, h_c) h_max = max(h_max, h_c) T_min = min(T_min, T_c) T_max = max(T_max, T_c) 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 face_flux = fields%face_flux(f) rho_c = transport%rho(owner) rho_p = rho_c if (allocated(fields%projection_rho)) then if (fields%projection_rho(owner) > 0.0_rk) rho_p = fields%projection_rho(owner) end if local_net_boundary_volume_flux = local_net_boundary_volume_flux + face_flux local_net_boundary_mass_flux_current_rho = local_net_boundary_mass_flux_current_rho + face_flux * rho_c local_net_boundary_mass_flux_projection_rho = local_net_boundary_mass_flux_projection_rho + face_flux * rho_p if (allocated(fields%mass_flux)) then mass_flux_stored = fields%mass_flux(f) local_net_boundary_mass_flux_stored = local_net_boundary_mass_flux_stored + mass_flux_stored end if end if end if end do end if sum_send = 0.0_rk sum_send(1) = local_volume sum_send(2) = local_mass_current sum_send(3) = local_mass_projection sum_send(4) = local_net_boundary_volume_flux sum_send(5) = local_net_boundary_mass_flux_stored sum_send(6) = local_net_boundary_mass_flux_current_rho sum_send(7) = local_net_boundary_mass_flux_projection_rho sum_send(8) = local_rho_current_l2_num sum_send(9) = local_rho_projection_l2_num sum_send(10) = local_rho_diff_l2_num sum_send(11) = local_rho_h_current sum_send(12) = local_rho_h_projection call MPI_Allreduce(sum_send, sum_recv, 12, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing transport conservation sums') min_send(1) = rho_current_min min_send(2) = rho_projection_min min_send(3) = h_min min_send(4) = T_min call MPI_Allreduce(min_send, min_recv, 4, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing transport conservation minima') max_send(1) = rho_current_max max_send(2) = rho_projection_max max_send(3) = rho_diff_max max_send(4) = h_max max_send(5) = T_max call MPI_Allreduce(max_send, max_recv, 5, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing transport conservation maxima') volume_total = sum_recv(1) mass_current = sum_recv(2) mass_projection = sum_recv(3) current_minus_projection_mass = mass_current - mass_projection net_boundary_volume_flux = sum_recv(4) net_boundary_mass_flux_stored = sum_recv(5) net_boundary_mass_flux_current_rho = sum_recv(6) net_boundary_mass_flux_projection_rho = sum_recv(7) stored_minus_current_rho_boundary_mass_flux = net_boundary_mass_flux_stored - net_boundary_mass_flux_current_rho stored_minus_projection_rho_boundary_mass_flux = net_boundary_mass_flux_stored - net_boundary_mass_flux_projection_rho rho_h_current = sum_recv(11) rho_h_projection = sum_recv(12) if (volume_total > 0.0_rk) then rho_current_mean = mass_current / volume_total rho_projection_mean = mass_projection / volume_total rho_current_l2 = sqrt(sum_recv(8) / volume_total) rho_projection_l2 = sqrt(sum_recv(9) / volume_total) rho_diff_l2 = sqrt(sum_recv(10) / volume_total) else rho_current_mean = 0.0_rk rho_projection_mean = 0.0_rk rho_current_l2 = 0.0_rk rho_projection_l2 = 0.0_rk rho_diff_l2 = 0.0_rk end if rho_current_min = min_recv(1) rho_projection_min = min_recv(2) h_min = min_recv(3) T_min = min_recv(4) rho_current_max = max_recv(1) rho_projection_max = max_recv(2) rho_diff_max = max_recv(3) h_max = max_recv(4) T_max = max_recv(5) if (transport_diag_initialized) then dt_prev = time - prev_time else dt_prev = 0.0_rk end if if (dt_prev > tiny(1.0_rk)) then delta_mass_current = mass_current - prev_mass_current mass_rate_current = delta_mass_current / dt_prev delta_mass_projection = mass_projection - prev_mass_projection mass_rate_projection = delta_mass_projection / dt_prev delta_rho_h_current = rho_h_current - prev_rho_h_current rho_h_rate_current = delta_rho_h_current / dt_prev else delta_mass_current = 0.0_rk mass_rate_current = 0.0_rk delta_mass_projection = 0.0_rk mass_rate_projection = 0.0_rk delta_rho_h_current = 0.0_rk rho_h_rate_current = 0.0_rk end if defect_current_stored = mass_rate_current + net_boundary_mass_flux_stored defect_current_current_rho = mass_rate_current + net_boundary_mass_flux_current_rho defect_projection_projection_rho = mass_rate_projection + net_boundary_mass_flux_projection_rho denom = max(max(abs(mass_rate_current), abs(net_boundary_mass_flux_stored)), tiny(1.0_rk)) rel_defect_current_stored = abs(defect_current_stored) / denom denom = max(max(abs(mass_rate_current), abs(net_boundary_mass_flux_current_rho)), tiny(1.0_rk)) rel_defect_current_current_rho = abs(defect_current_current_rho) / denom denom = max(max(abs(mass_rate_projection), abs(net_boundary_mass_flux_projection_rho)), tiny(1.0_rk)) rel_defect_projection_projection_rho = abs(defect_projection_projection_rho) / denom if (flow%rank == 0) then filename = trim(params%output_dir) // '/diagnostics/variable_density_transport_conservation.csv' if (.not. transport_diag_initialized) then open(newunit=unit_id, file=trim(filename), status='replace', action='write') write(unit_id,'(a)') 'step,time,volume_total,' // & 'mass_integral_current,mass_integral_projection,current_minus_projection_mass_integral,' // & 'net_boundary_volume_flux,net_boundary_mass_flux_stored,' // & 'net_boundary_mass_flux_current_rho,net_boundary_mass_flux_projection_rho,' // & 'stored_minus_current_rho_boundary_mass_flux,stored_minus_projection_rho_boundary_mass_flux,' // & 'delta_time_since_previous,delta_mass_current_since_previous,mass_rate_current_since_previous,' // & 'mass_balance_defect_current_stored_flux,relative_mass_balance_defect_current_stored_flux,' // & 'mass_balance_defect_current_current_rho_flux,relative_mass_balance_defect_current_current_rho_flux,' // & 'delta_mass_projection_since_previous,mass_rate_projection_since_previous,' // & 'mass_balance_defect_projection_projection_rho_flux,' // & 'relative_mass_balance_defect_projection_projection_rho_flux,' // & 'rho_current_min,rho_current_max,rho_current_mean,rho_current_l2,' // & 'rho_projection_min,rho_projection_max,rho_projection_mean,rho_projection_l2,' // & 'rho_current_minus_projection_max,rho_current_minus_projection_l2,' // & 'h_min,h_max,T_min,T_max,' // & 'rho_h_integral_current,rho_h_integral_projection,' // & 'delta_rho_h_current_since_previous,rho_h_current_rate_since_previous' else open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write') end if write(unit_id,'(i0,40(",",ES26.16E4))') step, time, volume_total, & mass_current, mass_projection, current_minus_projection_mass, & net_boundary_volume_flux, net_boundary_mass_flux_stored, & net_boundary_mass_flux_current_rho, net_boundary_mass_flux_projection_rho, & stored_minus_current_rho_boundary_mass_flux, stored_minus_projection_rho_boundary_mass_flux, & dt_prev, delta_mass_current, mass_rate_current, & defect_current_stored, rel_defect_current_stored, & defect_current_current_rho, rel_defect_current_current_rho, & delta_mass_projection, mass_rate_projection, & defect_projection_projection_rho, rel_defect_projection_projection_rho, & rho_current_min, rho_current_max, rho_current_mean, rho_current_l2, & rho_projection_min, rho_projection_max, rho_projection_mean, rho_projection_l2, & rho_diff_max, rho_diff_l2, h_min, h_max, T_min, T_max, & rho_h_current, rho_h_projection, delta_rho_h_current, rho_h_rate_current close(unit_id) end if transport_diag_initialized = .true. prev_time = time prev_mass_current = mass_current prev_mass_projection = mass_projection prev_rho_h_current = rho_h_current end subroutine write_variable_density_transport_conservation_diagnostics