Write conservative continuity residual diagnostics for variable-density mode.
This diagnostics-only routine separates the local density-history source relation from conservative mass continuity:
d(rho)/dt + div(rho u) = 0
and its expanded finite-volume form:
d(rho)/dt + rho div(u) + u.grad(rho) = 0
with u.grad(rho) estimated as div(rho u) - rho div(u).
| 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(transport_properties_t), | intent(in) | :: | transport | |||
| integer, | intent(in) | :: | step | |||
| real(kind=rk), | intent(in) | :: | time |
subroutine write_variable_density_continuity_residual_diagnostics(mesh, flow, params, fields, 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(transport_properties_t), intent(in) :: transport integer, intent(in) :: step real(rk), intent(in) :: time integer :: c, lf, fid, ierr, unit_id logical, save :: continuity_diag_initialized = .false. character(len=1024) :: filename real(rk) :: vol, rho_current, rho_projection real(rk) :: flux_out, mass_flux_out real(rk) :: divu, div_mass_flux, drho_dt real(rk) :: rho_divu, udotgradrho real(rk) :: history_res, conservative_res, expanded_res real(rk) :: local_volume, volume_total real(rk) :: local_integral_drho_dt_dV, integral_drho_dt_dV real(rk) :: local_integral_rho_divu_dV, integral_rho_divu_dV real(rk) :: local_integral_udotgradrho_dV, integral_udotgradrho_dV real(rk) :: local_integral_div_mass_flux_dV, integral_div_mass_flux_dV real(rk) :: local_integral_history_res_dV, integral_history_res_dV real(rk) :: local_integral_conservative_res_dV, integral_conservative_res_dV real(rk) :: local_integral_expanded_res_dV, integral_expanded_res_dV real(rk) :: local_history_res_max, history_res_max real(rk) :: local_conservative_res_max, conservative_res_max real(rk) :: local_expanded_res_max, expanded_res_max real(rk) :: local_udotgradrho_max, udotgradrho_max real(rk) :: local_history_res_l2_num, history_res_l2_num, history_res_l2 real(rk) :: local_conservative_res_l2_num, conservative_res_l2_num, conservative_res_l2 real(rk) :: local_expanded_res_l2_num, expanded_res_l2_num, expanded_res_l2 real(rk) :: local_udotgradrho_l2_num, udotgradrho_l2_num, udotgradrho_l2 real(rk) :: local_drho_dt_l2_num, drho_dt_l2_num, drho_dt_l2 real(rk) :: local_div_mass_flux_l2_num, div_mass_flux_l2_num, div_mass_flux_l2 real(rk) :: local_rho_divu_l2_num, rho_divu_l2_num, rho_divu_l2 real(rk) :: denom, rel_conservative_res_l2, rel_conservative_res_max real(rk) :: mpi_reduce_send(1), mpi_reduce_recv(1) if (.not. params%enable_variable_density) return if (.not. params%write_diagnostics) return if (.not. allocated(transport%rho)) return if (.not. allocated(fields%face_flux)) return local_volume = 0.0_rk local_integral_drho_dt_dV = 0.0_rk local_integral_rho_divu_dV = 0.0_rk local_integral_udotgradrho_dV = 0.0_rk local_integral_div_mass_flux_dV = 0.0_rk local_integral_history_res_dV = 0.0_rk local_integral_conservative_res_dV = 0.0_rk local_integral_expanded_res_dV = 0.0_rk local_history_res_max = 0.0_rk local_conservative_res_max = 0.0_rk local_expanded_res_max = 0.0_rk local_udotgradrho_max = 0.0_rk local_history_res_l2_num = 0.0_rk local_conservative_res_l2_num = 0.0_rk local_expanded_res_l2_num = 0.0_rk local_udotgradrho_l2_num = 0.0_rk local_drho_dt_l2_num = 0.0_rk local_div_mass_flux_l2_num = 0.0_rk local_rho_divu_l2_num = 0.0_rk do c = flow%first_cell, flow%last_cell vol = mesh%cells(c)%volume rho_current = transport%rho(c) rho_projection = rho_current if (allocated(fields%projection_rho)) rho_projection = fields%projection_rho(c) divu = 0.0_rk div_mass_flux = 0.0_rk 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) if (allocated(fields%mass_flux)) then mass_flux_out = fields%mass_flux(fid) else mass_flux_out = rho_projection * fields%face_flux(fid) end if else flux_out = -fields%face_flux(fid) if (allocated(fields%mass_flux)) then mass_flux_out = -fields%mass_flux(fid) else mass_flux_out = -rho_projection * fields%face_flux(fid) end if end if divu = divu + flux_out / vol div_mass_flux = div_mass_flux + mass_flux_out / vol end do if (params%dt > 0.0_rk) then drho_dt = (rho_current - rho_projection) / params%dt else drho_dt = 0.0_rk end if rho_divu = rho_current * divu udotgradrho = div_mass_flux - rho_divu history_res = drho_dt + rho_divu conservative_res = drho_dt + div_mass_flux expanded_res = drho_dt + rho_divu + udotgradrho local_volume = local_volume + vol local_integral_drho_dt_dV = local_integral_drho_dt_dV + drho_dt * vol local_integral_rho_divu_dV = local_integral_rho_divu_dV + rho_divu * vol local_integral_udotgradrho_dV = local_integral_udotgradrho_dV + udotgradrho * vol local_integral_div_mass_flux_dV = local_integral_div_mass_flux_dV + div_mass_flux * vol local_integral_history_res_dV = local_integral_history_res_dV + history_res * vol local_integral_conservative_res_dV = local_integral_conservative_res_dV + conservative_res * vol local_integral_expanded_res_dV = local_integral_expanded_res_dV + expanded_res * vol local_history_res_max = max(local_history_res_max, abs(history_res)) local_conservative_res_max = max(local_conservative_res_max, abs(conservative_res)) local_expanded_res_max = max(local_expanded_res_max, abs(expanded_res)) local_udotgradrho_max = max(local_udotgradrho_max, abs(udotgradrho)) local_history_res_l2_num = local_history_res_l2_num + history_res * history_res * vol local_conservative_res_l2_num = local_conservative_res_l2_num + conservative_res * conservative_res * vol local_expanded_res_l2_num = local_expanded_res_l2_num + expanded_res * expanded_res * vol local_udotgradrho_l2_num = local_udotgradrho_l2_num + udotgradrho * udotgradrho * vol local_drho_dt_l2_num = local_drho_dt_l2_num + drho_dt * drho_dt * vol local_div_mass_flux_l2_num = local_div_mass_flux_l2_num + div_mass_flux * div_mass_flux * vol local_rho_divu_l2_num = local_rho_divu_l2_num + rho_divu * rho_divu * vol end do 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 continuity volume') mpi_reduce_send(1) = local_integral_drho_dt_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_drho_dt_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity drho_dt') mpi_reduce_send(1) = local_integral_rho_divu_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_rho_divu_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity rho_divu') mpi_reduce_send(1) = local_integral_udotgradrho_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_udotgradrho_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity udotgradrho') mpi_reduce_send(1) = local_integral_div_mass_flux_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_div_mass_flux_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity div_mass_flux') mpi_reduce_send(1) = local_integral_history_res_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_history_res_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity history residual') mpi_reduce_send(1) = local_integral_conservative_res_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_conservative_res_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity conservative residual') mpi_reduce_send(1) = local_integral_expanded_res_dV call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) integral_expanded_res_dV = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity expanded residual') mpi_reduce_send(1) = local_history_res_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) history_res_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity history max') mpi_reduce_send(1) = local_conservative_res_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) conservative_res_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity conservative max') mpi_reduce_send(1) = local_expanded_res_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) expanded_res_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity expanded max') mpi_reduce_send(1) = local_udotgradrho_max call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr) udotgradrho_max = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity udotgradrho max') mpi_reduce_send(1) = local_history_res_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) history_res_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity history l2') mpi_reduce_send(1) = local_conservative_res_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) conservative_res_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity conservative l2') mpi_reduce_send(1) = local_expanded_res_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) expanded_res_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity expanded l2') mpi_reduce_send(1) = local_udotgradrho_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) udotgradrho_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity udotgradrho l2') mpi_reduce_send(1) = local_drho_dt_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) drho_dt_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity drho_dt l2') mpi_reduce_send(1) = local_div_mass_flux_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) div_mass_flux_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity div_mass_flux l2') mpi_reduce_send(1) = local_rho_divu_l2_num call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) rho_divu_l2_num = mpi_reduce_recv(1) if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing continuity rho_divu l2') if (volume_total > 0.0_rk) then history_res_l2 = sqrt(history_res_l2_num / volume_total) conservative_res_l2 = sqrt(conservative_res_l2_num / volume_total) expanded_res_l2 = sqrt(expanded_res_l2_num / volume_total) udotgradrho_l2 = sqrt(udotgradrho_l2_num / volume_total) drho_dt_l2 = sqrt(drho_dt_l2_num / volume_total) div_mass_flux_l2 = sqrt(div_mass_flux_l2_num / volume_total) rho_divu_l2 = sqrt(rho_divu_l2_num / volume_total) else history_res_l2 = 0.0_rk conservative_res_l2 = 0.0_rk expanded_res_l2 = 0.0_rk udotgradrho_l2 = 0.0_rk drho_dt_l2 = 0.0_rk div_mass_flux_l2 = 0.0_rk rho_divu_l2 = 0.0_rk end if denom = max(div_mass_flux_l2 + drho_dt_l2, tiny(1.0_rk)) rel_conservative_res_l2 = conservative_res_l2 / denom denom = max(max(abs(integral_div_mass_flux_dV), abs(integral_drho_dt_dV)), tiny(1.0_rk)) rel_conservative_res_max = abs(integral_conservative_res_dV) / denom if (flow%rank == 0) then filename = trim(params%output_dir) // '/diagnostics/variable_density_continuity_residual.csv' if (.not. continuity_diag_initialized) then open(newunit=unit_id, file=trim(filename), status='replace', action='write') continuity_diag_initialized = .true. write(unit_id,'(a)') 'step,time,volume_total,' // & 'integral_drho_dt_dV,integral_rho_current_divu_dV,' // & 'integral_udotgradrho_dV,integral_div_mass_flux_dV,' // & 'integral_drho_dt_plus_rho_current_divu_dV,' // & 'integral_drho_dt_plus_div_mass_flux_dV,' // & 'integral_drho_dt_plus_rho_current_divu_plus_udotgradrho_dV,' // & 'history_residual_max,history_residual_l2,' // & 'conservative_residual_max,conservative_residual_l2,' // & 'expanded_residual_max,expanded_residual_l2,' // & 'udotgradrho_max,udotgradrho_l2,' // & 'drho_dt_l2,div_mass_flux_l2,rho_current_divu_l2,' // & 'relative_conservative_residual_l2,' // & 'relative_integral_conservative_residual' else open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write') end if write(unit_id,'(i0,22(",",ES26.16E4))') step, time, volume_total, & integral_drho_dt_dV, integral_rho_divu_dV, integral_udotgradrho_dV, & integral_div_mass_flux_dV, integral_history_res_dV, & integral_conservative_res_dV, integral_expanded_res_dV, & history_res_max, history_res_l2, conservative_res_max, conservative_res_l2, & expanded_res_max, expanded_res_l2, udotgradrho_max, udotgradrho_l2, & drho_dt_l2, div_mass_flux_l2, rho_divu_l2, rel_conservative_res_l2, & rel_conservative_res_max close(unit_id) end if end subroutine write_variable_density_continuity_residual_diagnostics