diff --git a/CMakeLists.txt b/CMakeLists.txt index 0dca30c..cc9519e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ # # This file is part of Theseus. # -# SPDX-License-Identifier: MIT +# SPDX-License-Identifier: BSD-3-Clause cmake_minimum_required(VERSION 3.15) set(CMAKE_CXX_STANDARD 17) @@ -22,6 +22,10 @@ option(ENABLE_TIMER_OUTPUT_ALLRANKS "Output timers on all ranks." OFF) option(ENABLE_TIMER_SYNC_DEVICE "Fence timers with device sync." OFF) option(ENABLE_CUDA "Compile CUDA kernels" OFF) option(THESEUS_WITH_PLATO "Enable PLATO API for thermo tables" OFF) +option(ENABLE_POINT_PARALLEL_INTERIOR_FACES + "Use one device thread per interior-face point." OFF) +option(ENABLE_POINT_PARALLEL_VOLUME + "Use one device thread per element volume point where supported." OFF) ## NOTE: PLATO_DB_PATH can be specified at build-time to customize ## the location of the database used by PLATO. Tries to default ## to a reasonable location. @@ -142,20 +146,32 @@ else() endif() if(ENABLE_TIMERS) - message(STATUS "Enabling performance timers.") + message(STATUS "Enabling detailed performance timers.") target_compile_definitions(theseus PRIVATE ENABLE_TIMERS) - if(ENABLE_TIMER_BARRIER) - message(STATUS " - Timing with barriers") - target_compile_definitions(theseus PRIVATE TIMER_BARRIER) - endif() if(ENABLE_TIMER_OUTPUT_ALLRANKS) message(STATUS " - Timing output on all ranks") target_compile_definitions(theseus PRIVATE TIMER_OUTPUT_ALLRANKS) endif() - if(ENABLE_TIMER_SYNC_DEVICE) - message(STATUS " - Timing device sync enabled") - target_compile_definitions(theseus PRIVATE TIMER_SYNC_DEVICE) - endif() +else() + message(STATUS "Detailed performance timers disabled.") +endif() +if(ENABLE_TIMER_SYNC_DEVICE) + message(STATUS " - Timing device sync enabled") + target_compile_definitions(theseus PRIVATE TIMER_SYNC_DEVICE) +endif() +if(ENABLE_TIMER_BARRIER) + message(STATUS " - Timing barriers enabled") + target_compile_definitions(theseus PRIVATE TIMER_BARRIER) +endif() + +if(ENABLE_POINT_PARALLEL_INTERIOR_FACES) + message(STATUS "Using point-parallel interior-face assembly") + target_compile_definitions(theseus PRIVATE POINT_PARALLEL_INTERIOR_FACES) +endif() + +if(ENABLE_POINT_PARALLEL_VOLUME) + message(STATUS "Using point-parallel volume assembly where supported") + target_compile_definitions(theseus PRIVATE POINT_PARALLEL_VOLUME) endif() if(SUBCELL_FV_BLENDING) diff --git a/include/DGSEMIntegrator.hpp b/include/DGSEMIntegrator.hpp index f2751ef..83b4b51 100644 --- a/include/DGSEMIntegrator.hpp +++ b/include/DGSEMIntegrator.hpp @@ -123,41 +123,212 @@ namespace Theseus return max_char_speed; } + template + MFEM_HOST_DEVICE inline + static mfem::real_t AssembleVolumePointKernel( + const ContextType &ctx, const mfem::real_t *el_u, + const mfem::real_t *elJac_d, const mfem::real_t *elMetric_d, + const int point, mfem::real_t *el_dudt) + { + const int Np_x = ctx.Np_x; + const int Np_y = ctx.Np_y; + const int Np_z = ctx.Np_z; + const int dim = ctx.dim; + const int neq = ctx.num_equations; + const int dof = ctx.ndof_scalar_el; + const mfem::real_t *Dhat2_d = ctx.Dhat2_d; + const int i = point % Np_x; + const int j = (point / Np_x) % Np_y; + const int k = point / (Np_x*Np_y); + + mfem::real_t state_lower[Theseus::MAXEQ]; + mfem::real_t state_upper[Theseus::MAXEQ]; + mfem::real_t flux[Theseus::MAXEQ]; + mfem::real_t point_rate[Theseus::MAXEQ] = {0.0}; + mfem::real_t max_char_speed = 0.0; + + for (int m = 0; m < Np_x; ++m) + { + if (m == i) { continue; } + const int lower = m < i ? m : i; + const int upper = m < i ? i : m; + const int lower_point = k*Np_y*Np_x + j*Np_x + lower; + const int upper_point = k*Np_y*Np_x + j*Np_x + upper; + Kernels::el_gather_state(el_u, dof, neq, lower_point, state_lower); + Kernels::el_gather_state(el_u, dof, neq, upper_point, state_upper); + const mfem::real_t char_speed = ctx.iflux.ComputeVolumeFlux( + ctx.gas, state_lower, state_upper, + elMetric_d + lower_point*dim*dim, + elMetric_d + upper_point*dim*dim, flux); + max_char_speed = Kernels::rmax(max_char_speed, char_speed); + const mfem::real_t coefficient = Dhat2_d[m + Np_x*i]; + for (int q = 0; q < neq; ++q) + { + point_rate[q] += coefficient*flux[q]; + } + } + + if (dim > 1) + { + for (int m = 0; m < Np_y; ++m) + { + if (m == j) { continue; } + const int lower = m < j ? m : j; + const int upper = m < j ? j : m; + const int lower_point = k*Np_y*Np_x + lower*Np_x + i; + const int upper_point = k*Np_y*Np_x + upper*Np_x + i; + Kernels::el_gather_state( + el_u, dof, neq, lower_point, state_lower); + Kernels::el_gather_state( + el_u, dof, neq, upper_point, state_upper); + const mfem::real_t char_speed = ctx.iflux.ComputeVolumeFlux( + ctx.gas, state_lower, state_upper, + elMetric_d + lower_point*dim*dim + dim, + elMetric_d + upper_point*dim*dim + dim, flux); + max_char_speed = Kernels::rmax(max_char_speed, char_speed); + const mfem::real_t coefficient = Dhat2_d[m + Np_y*j]; + for (int q = 0; q < neq; ++q) + { + point_rate[q] += coefficient*flux[q]; + } + } + } + + if (dim > 2) + { + for (int m = 0; m < Np_z; ++m) + { + if (m == k) { continue; } + const int lower = m < k ? m : k; + const int upper = m < k ? k : m; + const int lower_point = lower*Np_y*Np_x + j*Np_x + i; + const int upper_point = upper*Np_y*Np_x + j*Np_x + i; + Kernels::el_gather_state( + el_u, dof, neq, lower_point, state_lower); + Kernels::el_gather_state( + el_u, dof, neq, upper_point, state_upper); + const mfem::real_t char_speed = ctx.iflux.ComputeVolumeFlux( + ctx.gas, state_lower, state_upper, + elMetric_d + lower_point*dim*dim + 2*dim, + elMetric_d + upper_point*dim*dim + 2*dim, flux); + max_char_speed = Kernels::rmax(max_char_speed, char_speed); + const mfem::real_t coefficient = Dhat2_d[m + Np_z*k]; + for (int q = 0; q < neq; ++q) + { + point_rate[q] += coefficient*flux[q]; + } + } + } + + Kernels::el_scatter_assign( + point_rate, dof, neq, point, -1.0/elJac_d[point], el_dudt); + return max_char_speed; + } + + template + MFEM_HOST_DEVICE static mfem::real_t AssembleFacePointKernel(const ContextT &ctx, + const mfem::real_t *u_face, + const mfem::real_t *nor_point, + const mfem::real_t w_minus, + const mfem::real_t w_plus, + const int fp, + mfem::real_t *rhs_face) + { + mfem::real_t point_flux[Theseus::MAXEQ]; + mfem::real_t qMinus[Theseus::MAXEQ]; + mfem::real_t qPlus[Theseus::MAXEQ]; + const int neq = ctx.num_equations; + for(int q = 0; q < neq; ++q){ + qMinus[q] = u_face[ctx.iface_idx(0, fp, q)]; + qPlus[q] = u_face[ctx.iface_idx(1, fp, q)]; + } + + const mfem::real_t char_speed = + ctx.iflux.ComputeFaceFlux(ctx.gas, qMinus, qPlus, nor_point, point_flux); + + for(int q = 0; q < neq; ++q){ + rhs_face[ctx.iface_idx(0, fp, q)] = -w_minus * point_flux[q]; + rhs_face[ctx.iface_idx(1, fp, q)] = w_plus * point_flux[q]; + } + + return char_speed; + } + template MFEM_HOST_DEVICE static mfem::real_t AssembleElementFaceKernel(const ContextT &ctx, const mfem::real_t *u_face, const mfem::real_t *nor_face,const mfem::real_t *w_minus, const mfem::real_t *w_plus, mfem::real_t *rhs_face) { mfem::real_t max_char_speed = 0.0; + const int nfp = ctx.num_face_points; + const int dim = ctx.dim; + for (int fp = 0; fp < nfp; ++fp) + { + const mfem::real_t char_speed = + AssembleFacePointKernel(ctx, u_face, nor_face + fp*dim, + w_minus[fp], w_plus[fp], fp, rhs_face); + max_char_speed = Kernels::rmax(max_char_speed, char_speed); + } + return max_char_speed; + } + + template + MFEM_HOST_DEVICE static mfem::real_t AssembleViscousFacePointKernel( + const ContextT &ctx, const mfem::real_t *u_face, + const mfem::real_t *nor_point, const mfem::real_t w_minus, + const mfem::real_t w_plus, const mfem::real_t *dprim_face_x, + const mfem::real_t *dprim_face_y, const mfem::real_t *dprim_face_z, + const mfem::real_t radius, const int fp, mfem::real_t *rhs_face) + { mfem::real_t point_flux[Theseus::MAXEQ]; + mfem::real_t vflux_minus[Theseus::MAXEQ][Theseus::MAXDIM]; + mfem::real_t vflux_plus[Theseus::MAXEQ][Theseus::MAXDIM]; mfem::real_t qMinus[Theseus::MAXEQ]; mfem::real_t qPlus[Theseus::MAXEQ]; - const int nfp = ctx.num_face_points; + mfem::real_t gradPrim_plus[Theseus::MAXDIM][Theseus::MAXEQ]; + mfem::real_t gradPrim_minus[Theseus::MAXDIM][Theseus::MAXEQ]; + const mfem::real_t *dprim_face[Theseus::MAXDIM] = { + dprim_face_x, dprim_face_y, dprim_face_z}; const int neq = ctx.num_equations; const int dim = ctx.dim; - // auto idx = [=](int side, int fp, int eq) -> int - // { - // return (((side)*neq + eq)*nfp + fp); - // }; - for (int i = 0; i < nfp; i++) - { - const mfem::real_t *nor_d = nor_face + i*dim; - const mfem::real_t wminus = -w_minus[i]; - const mfem::real_t wplus = w_plus[i]; - // Could avoid these copy-in,out - for(int j = 0;j < neq;j++){ - qMinus[j] = u_face[ctx.iface_idx(0,i,j)]; - qPlus[j] = u_face[ctx.iface_idx(1,i,j)]; - } - max_char_speed = \ - Kernels::rmax(max_char_speed, ctx.iflux.ComputeFaceFlux(ctx.gas, qMinus, qPlus, - nor_d, point_flux)); - for(int j = 0;j < neq;j++){ - rhs_face[ctx.iface_idx(0, i, j)] = wminus * point_flux[j]; - rhs_face[ctx.iface_idx(1, i, j)] = wplus * point_flux[j]; - } + + for(int q = 0; q < neq; ++q){ + const int minus_index = ctx.iface_idx(0, fp, q); + const int plus_index = ctx.iface_idx(1, fp, q); + qMinus[q] = u_face[minus_index]; + qPlus[q] = u_face[plus_index]; + for(int idim = 0; idim < dim; ++idim){ + gradPrim_minus[idim][q] = dprim_face[idim][minus_index]; + gradPrim_plus[idim][q] = dprim_face[idim][plus_index]; } - return max_char_speed; + } + + const mfem::real_t char_speed = + ctx.iflux.ComputeFaceFlux(ctx.gas, qMinus, qPlus, nor_point, point_flux); + + NavierStokesFlux::ComputeViscousFluxKernel( + ctx.gas, qMinus, gradPrim_minus[0], gradPrim_minus[1], + gradPrim_minus[2], vflux_minus, ctx.axisymmetric, + ctx.axisymmetric ? radius : 0.0); + NavierStokesFlux::ComputeViscousFluxKernel( + ctx.gas, qPlus, gradPrim_plus[0], gradPrim_plus[1], + gradPrim_plus[2], vflux_plus, ctx.axisymmetric, + ctx.axisymmetric ? radius : 0.0); + + for(int q = 0; q < neq; ++q){ + for(int idim = 0; idim < dim; ++idim){ + const mfem::real_t avg = + 0.5*(vflux_minus[q][idim] + vflux_plus[q][idim]); + point_flux[q] -= nor_point[idim]*avg; + } + } + + for(int q = 0; q < neq; ++q){ + rhs_face[ctx.iface_idx(0, fp, q)] = -w_minus * point_flux[q]; + rhs_face[ctx.iface_idx(1, fp, q)] = w_plus * point_flux[q]; + } + + return char_speed; } template @@ -169,76 +340,16 @@ namespace Theseus mfem::real_t *rhs_face) { mfem::real_t max_char_speed = 0.0; - mfem::real_t point_flux[Theseus::MAXEQ]; - mfem::real_t vflux_minus[Theseus::MAXEQ][Theseus::MAXDIM]; - mfem::real_t vflux_plus[Theseus::MAXEQ][Theseus::MAXDIM]; - mfem::real_t qMinus[Theseus::MAXEQ]; - mfem::real_t qPlus[Theseus::MAXEQ]; - mfem::real_t gradPrim_plus[Theseus::MAXDIM][Theseus::MAXEQ]; - mfem::real_t gradPrim_minus[Theseus::MAXDIM][Theseus::MAXEQ]; - const mfem::real_t *dprim_face[Theseus::MAXDIM] = {dprim_face_x, dprim_face_y, dprim_face_z}; const int nfp = ctx.num_face_points; - const int neq = ctx.num_equations; const int dim = ctx.dim; - // auto idx = [=](int side, int fp, int eq) -> int - // { - // return (((side)*neq + eq)*nfp + fp); - // }; - for (int i = 0; i < nfp; i++) + for (int fp = 0; fp < nfp; ++fp) { - const mfem::real_t *nor_d = nor_face + i*dim; - const mfem::real_t wminus = -w_minus[i]; - const mfem::real_t wplus = w_plus[i]; - // Could avoid these copy-in,out - for(int j = 0;j < neq;j++){ - int minus_index = ctx.iface_idx(0, i, j); - int plus_index = ctx.iface_idx(1, i, j); - qMinus[j] = u_face[minus_index]; - qPlus[j] = u_face[plus_index]; - for(int idim = 0;idim < dim;idim++){ - gradPrim_minus[idim][j] = dprim_face[idim][minus_index]; - gradPrim_plus[idim][j] = dprim_face[idim][plus_index]; - } - } - max_char_speed = \ - Kernels::rmax(max_char_speed, ctx.iflux.ComputeFaceFlux(ctx.gas, qMinus, qPlus, - nor_d, point_flux)); - - // Here, point_flux is +(F_inv * Normal) - - // Grab the viscous flux - NavierStokesFlux::ComputeViscousFluxKernel(ctx.gas, qMinus, - gradPrim_minus[0], - gradPrim_minus[1], - gradPrim_minus[2], vflux_minus, - ctx.axisymmetric, - ctx.axisymmetric ? face_radius[i] : 0.0); - NavierStokesFlux::ComputeViscousFluxKernel(ctx.gas, qPlus, - gradPrim_plus[0], - gradPrim_plus[1], - gradPrim_plus[2], vflux_plus, - ctx.axisymmetric, - ctx.axisymmetric ? face_radius[i] : 0.0); - - // Now we have vflux(+) and vflux(-) - // In this loop: - // - average vflux - // - dot avg vflux with nor - // - accumulate dotted (avg*n) into point_flux - for(int j = 0;j < neq;j++){ - for(int idim = 0;idim < dim;idim++){ - mfem::real_t avg = 0.5*(vflux_minus[j][idim] + vflux_plus[j][idim]); - point_flux[j] -= nor_d[idim]*avg; - } - } - // So now: point_flux = +(F_inv * Normal) -(F^bar_visc * Normal) - // in this loop: - // - SET/Overwrite rhs_face - // - NEGATE the (-) face point_flux to properly orient - for(int j = 0;j < neq;j++){ - rhs_face[ctx.iface_idx(0, i, j)] = wminus * point_flux[j]; - rhs_face[ctx.iface_idx(1, i, j)] = wplus * point_flux[j]; - } + const mfem::real_t radius = + ctx.axisymmetric ? face_radius[fp] : 0.0; + const mfem::real_t char_speed = AssembleViscousFacePointKernel( + ctx, u_face, nor_face + fp*dim, w_minus[fp], w_plus[fp], + dprim_face_x, dprim_face_y, dprim_face_z, radius, fp, rhs_face); + max_char_speed = Kernels::rmax(max_char_speed, char_speed); } return max_char_speed; } @@ -401,6 +512,117 @@ namespace Theseus } + template + MFEM_HOST_DEVICE inline + static void AssembleViscousVolumePointKernel( + const ContextType &ctx, const mfem::real_t *el_u, + const mfem::real_t *elJac_d, const mfem::real_t *elMetric_d, + const mfem::real_t *elRadius_d, + const mfem::real_t *el_gradprim_x, + const mfem::real_t *el_gradprim_y, + const mfem::real_t *el_gradprim_z, + const int point, mfem::real_t *el_dudt) + { + const int Np_x = ctx.Np_x; + const int Np_y = ctx.Np_y; + const int Np_z = ctx.Np_z; + const int dim = ctx.dim; + const int neq = ctx.num_equations; + const int dof = ctx.ndof_scalar_el; + const mfem::real_t *Dhat_d = ctx.Dhat_d; + const int i = point % Np_x; + const int j = (point / Np_x) % Np_y; + const int k = point / (Np_x*Np_y); + + mfem::real_t state[Theseus::MAXEQ] = {0.0}; + mfem::real_t dqx[Theseus::MAXEQ] = {0.0}; + mfem::real_t dqy[Theseus::MAXEQ] = {0.0}; + mfem::real_t dqz[Theseus::MAXEQ] = {0.0}; + mfem::real_t f_ref[Theseus::MAXEQ] = {0.0}; + mfem::real_t dU_viscous[Theseus::MAXEQ] = {0.0}; + + for (int l = 0; l < Np_x; ++l) + { + const int sample = k*Np_y*Np_x + j*Np_x + l; + const mfem::real_t coefficient = Dhat_d[l + Np_x*i]; + Kernels::el_gather_state(el_u, dof, neq, sample, state); + Kernels::el_gather_grad_state( + el_gradprim_x, el_gradprim_y, el_gradprim_z, dim, dof, neq, + sample, dqx, dqy, dqz); + Theseus::NavierStokesFlux::compute_ref_viscous_flux( + ctx.gas, dim, neq, state, dqx, dqy, dqz, + elMetric_d + sample*dim*dim, f_ref, ctx.axisymmetric, + ctx.axisymmetric ? elRadius_d[sample] : 0.0); + for (int q = 0; q < neq; ++q) + { + dU_viscous[q] += coefficient*f_ref[q]; + } + } + + if (dim > 1) + { + for (int l = 0; l < Np_y; ++l) + { + const int sample = k*Np_y*Np_x + l*Np_x + i; + const mfem::real_t coefficient = Dhat_d[l + Np_y*j]; + Kernels::el_gather_state(el_u, dof, neq, sample, state); + Kernels::el_gather_grad_state( + el_gradprim_x, el_gradprim_y, el_gradprim_z, dim, dof, neq, + sample, dqx, dqy, dqz); + Theseus::NavierStokesFlux::compute_ref_viscous_flux( + ctx.gas, dim, neq, state, dqx, dqy, dqz, + elMetric_d + sample*dim*dim + dim, f_ref, + ctx.axisymmetric, + ctx.axisymmetric ? elRadius_d[sample] : 0.0); + for (int q = 0; q < neq; ++q) + { + dU_viscous[q] += coefficient*f_ref[q]; + } + } + } + + if (dim > 2) + { + for (int l = 0; l < Np_z; ++l) + { + const int sample = l*Np_y*Np_x + j*Np_x + i; + const mfem::real_t coefficient = Dhat_d[l + Np_z*k]; + Kernels::el_gather_state(el_u, dof, neq, sample, state); + Kernels::el_gather_grad_state( + el_gradprim_x, el_gradprim_y, el_gradprim_z, dim, dof, neq, + sample, dqx, dqy, dqz); + Theseus::NavierStokesFlux::compute_ref_viscous_flux( + ctx.gas, dim, neq, state, dqx, dqy, dqz, + elMetric_d + sample*dim*dim + 2*dim, f_ref, + ctx.axisymmetric, + ctx.axisymmetric ? elRadius_d[sample] : 0.0); + for (int q = 0; q < neq; ++q) + { + dU_viscous[q] += coefficient*f_ref[q]; + } + } + } + + Kernels::el_scatter_add( + dU_viscous, dof, neq, point, 1.0/elJac_d[point], el_dudt); + if (ctx.axisymmetric) + { + Kernels::el_gather_state(el_u, dof, neq, point, state); + Kernels::el_gather_grad_state( + el_gradprim_x, el_gradprim_y, el_gradprim_z, dim, dof, neq, + point, dqx, dqy, dqz); + mfem::real_t source[Theseus::MAXEQ] = {0.0}; + if (!AddAxisymmetricViscousSourceAwayFromAxis( + ctx.gas, state, dqx, dqy, dqz, elRadius_d[point], source)) + { + AddAxisymmetricViscousSourceAtAxis( + ctx, el_u, el_gradprim_x, el_gradprim_y, el_gradprim_z, + elRadius_d, elJac_d, elMetric_d, point, source); + } + Kernels::el_scatter_add(source, dof, neq, point, 1.0, el_dudt); + } + } + template MFEM_HOST_DEVICE inline static void AssembleViscousElementVolumeKernel(const ContextType &ctx, @@ -541,6 +763,94 @@ namespace Theseus } } + template + MFEM_HOST_DEVICE inline + static void AssembleGradVolumePointKernel( + const ContextType &ctx, const mfem::real_t *el_u, + const mfem::real_t *elJac_d, const mfem::real_t *elMetric_d, + const int point, mfem::real_t *el_grad_u[Theseus::MAXDIM]) + { + const int Np_x = ctx.Np_x; + const int Np_y = ctx.Np_y; + const int neq = ctx.num_equations; + const int dim = ctx.dim; + const int dof = ctx.ndof_scalar_el; + const mfem::real_t *D_d = ctx.D_d; + + const int i = point % Np_x; + const int j = (point / Np_x) % Np_y; + const int k = point / (Np_x * Np_y); + + mfem::real_t dudxi[Theseus::MAXEQ] = {0.0}; + mfem::real_t dudeta[Theseus::MAXEQ] = {0.0}; + mfem::real_t dudzeta[Theseus::MAXEQ] = {0.0}; + + for (int l = 0; l < Np_x; ++l) + { + const int sample = k*Np_y*Np_x + j*Np_x + l; + const mfem::real_t coefficient = D_d[l + Np_x*i]; + for (int q = 0; q < neq; ++q) + { + dudxi[q] += el_u[sample + q*dof] * coefficient; + } + } + + if (dim > 1) + { + for (int l = 0; l < Np_y; ++l) + { + const int sample = k*Np_y*Np_x + l*Np_x + i; + const mfem::real_t coefficient = D_d[l + Np_y*j]; + for (int q = 0; q < neq; ++q) + { + dudeta[q] += el_u[sample + q*dof] * coefficient; + } + } + } + + if (dim > 2) + { + for (int l = 0; l < ctx.Np_z; ++l) + { + const int sample = l*Np_y*Np_x + j*Np_x + i; + const mfem::real_t coefficient = D_d[l + ctx.Np_z*k]; + for (int q = 0; q < neq; ++q) + { + dudzeta[q] += el_u[sample + q*dof] * coefficient; + } + } + } + + const mfem::real_t invJ = 1.0 / elJac_d[point]; + const mfem::real_t *adj = elMetric_d + point*dim*dim; + for (int q = 0; q < neq; ++q) + { + if (dim == 1) + { + el_grad_u[0][point + q*dof] = invJ*dudxi[q]*adj[0]; + } + else if (dim == 2) + { + el_grad_u[0][point + q*dof] = + invJ*(dudxi[q]*adj[0] + dudeta[q]*adj[2]); + el_grad_u[1][point + q*dof] = + invJ*(dudxi[q]*adj[1] + dudeta[q]*adj[3]); + } + else + { + el_grad_u[0][point + q*dof] = + invJ*(dudxi[q]*adj[0] + dudeta[q]*adj[3] + + dudzeta[q]*adj[6]); + el_grad_u[1][point + q*dof] = + invJ*(dudxi[q]*adj[1] + dudeta[q]*adj[4] + + dudzeta[q]*adj[7]); + el_grad_u[2][point + q*dof] = + invJ*(dudxi[q]*adj[2] + dudeta[q]*adj[5] + + dudzeta[q]*adj[8]); + } + } + } + template MFEM_HOST_DEVICE inline static void AssembleGradElementVolumeKernel(const ContextType &ctx, @@ -702,45 +1012,56 @@ namespace Theseus template MFEM_HOST_DEVICE inline - static void AssembleGradInteriorFaceKernel(const ContextT &ctx, + static void AssembleGradInteriorFacePointKernel( + const ContextT &ctx, const mfem::real_t *u_face, - const mfem::real_t *nor_face, - const mfem::real_t *w_minus, - const mfem::real_t *w_plus, + const mfem::real_t *nor_point, + const mfem::real_t w_minus, + const mfem::real_t w_plus, + const int fp, mfem::real_t *rhs_face[Theseus::MAXDIM]) { - const int nfp = ctx.num_face_points; const int neq = ctx.num_equations; const int dim = ctx.dim; - mfem::real_t qMinus[Theseus::MAXEQ]; - mfem::real_t qPlus[Theseus::MAXEQ]; mfem::real_t jump[Theseus::MAXEQ]; - for (int i = 0; i < nfp; ++i) + for (int q = 0; q < neq; ++q) { - const mfem::real_t *nor_d = nor_face + i * dim; + jump[q] = mfem::real_t(0.5) * + (u_face[ctx.iface_idx(1, fp, q)] - + u_face[ctx.iface_idx(0, fp, q)]); + } - const mfem::real_t wminus = w_minus[i]; - const mfem::real_t wplus = w_plus[i]; + for (int idim = 0; idim < dim; ++idim){ + mfem::real_t *rhs_d = rhs_face[idim]; + const mfem::real_t n_d = nor_point[idim]; + for (int q = 0; q < neq; ++q) + { + const mfem::real_t f_d = jump[q]*n_d; + rhs_d[ctx.iface_idx(0, fp, q)] = w_minus * f_d; + rhs_d[ctx.iface_idx(1, fp, q)] = w_plus * f_d; + } + } + } - for (int q = 0; q < neq; ++q) - { - qMinus[q] = u_face[ctx.iface_idx(0, i, q)]; - qPlus[q] = u_face[ctx.iface_idx(1, i, q)]; - jump[q] = mfem::real_t(0.5) * (qPlus[q] - qMinus[q]); - } + template + MFEM_HOST_DEVICE inline + static void AssembleGradInteriorFaceKernel(const ContextT &ctx, + const mfem::real_t *u_face, + const mfem::real_t *nor_face, + const mfem::real_t *w_minus, + const mfem::real_t *w_plus, + mfem::real_t *rhs_face[Theseus::MAXDIM]) + { + const int nfp = ctx.num_face_points; + const int dim = ctx.dim; - for ( int idim = 0;idim < dim;idim++){ - mfem::real_t *rhs_d = rhs_face[idim]; - const mfem::real_t n_d = nor_d[idim]; - for (int q = 0; q < neq; ++q) - { - const mfem::real_t f_d = jump[q]*n_d; - rhs_d[ctx.iface_idx(0, i, q)] = wminus * f_d; - rhs_d[ctx.iface_idx(1, i, q)] = wplus * f_d; - } - } + for (int fp = 0; fp < nfp; ++fp) + { + AssembleGradInteriorFacePointKernel( + ctx, u_face, nor_face + fp*dim, w_minus[fp], w_plus[fp], fp, + rhs_face); } } diff --git a/include/EulerOperator_impl.hpp b/include/EulerOperator_impl.hpp index 49b4f84..22419fa 100644 --- a/include/EulerOperator_impl.hpp +++ b/include/EulerOperator_impl.hpp @@ -102,6 +102,68 @@ namespace Theseus mfem::real_t *ws_d = dc.elWaveSpeed_d; // Inside the FORALL below, executed on device +#ifdef POINT_PARALLEL_VOLUME + const int npoints = ne * ndof; + mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p) + { + const int e = p / ndof; + const int point = p % ndof; + const int attr = elem_attr_d[e]; + if (attr_marker_d[attr-1] == 0) { + ws_d[p] = 0.0; + return; + } + + const int element_offset = e * estride; + ws_d[p] = DGSEMIntegrator::AssembleVolumePointKernel( + dc, Ue_d + element_offset, elJac_d + e*jac_stride, + elMetric_d + e*metric_stride, point, dUe_d + element_offset); + }); + + mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e) + { + const int attr = elem_attr_d[e]; + if (attr_marker_d[attr-1] == 0) { + return; + } + + const int element_offset = e * estride; + const mfem::real_t *u_el = Ue_d + element_offset; + mfem::real_t *du_el = dUe_d + element_offset; + const mfem::real_t *jac_el = elJac_d + e*jac_stride; + const mfem::real_t *metric_el = elMetric_d + e*metric_stride; + const mfem::real_t *radius_el = + dc.axisymmetric ? elRadius_d + e*jac_stride : nullptr; + mfem::real_t char_speed = ws_d[e*ndof]; + +#ifdef SUBCELL_FV_BLENDING + const mfem::real_t alpha_fv = alpha_d[e]; + if (alpha_fv > 1e-16) { + const mfem::real_t alpha_dg = 1.0 - alpha_fv; + mfem::real_t *du_fv = dUfv_d + element_offset; + const mfem::real_t *el_metric_xi = + metric_xi_d + e*npe_metric_xi*dim; + const mfem::real_t *el_metric_eta = dim > 1 ? + metric_eta_d + e*npe_metric_eta*dim : nullptr; + const mfem::real_t *el_metric_zeta = dim > 2 ? + metric_zeta_d + e*npe_metric_zeta*dim : nullptr; + const mfem::real_t fv_char_speed = + DGSEMIntegrator::ComputeFVFluxesKernel( + dc, u_el, jac_el, el_metric_xi, el_metric_eta, + el_metric_zeta, du_fv); + for (int value = 0; value < estride; ++value) { + du_el[value] = + alpha_dg*du_el[value] + alpha_fv*du_fv[value]; + } + char_speed = Kernels::rmax(char_speed, fv_char_speed); + } +#endif + + AddAxisymmetricEulerElementSource( + dc, u_el, radius_el, jac_el, metric_el, du_el); + ws_d[e*ndof] = char_speed; + }); +#else mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e) { @@ -152,6 +214,7 @@ namespace Theseus ws_d[e] = cs_el; }); +#endif // Scatter RHS back to storage operator_cache.restr_v->AddMultTranspose(dUe, pdudt); @@ -160,9 +223,14 @@ namespace Theseus // - Reduce for rank-local max_char_speed const mfem::real_t *ws = operator_cache.elWaveSpeed.HostRead(); mfem::real_t max_char_speed = 0.0; - for(int e = 0;e < operator_cache.num_elements;e++) +#ifdef POINT_PARALLEL_VOLUME + const int num_wave_speeds = ne * ndof; +#else + const int num_wave_speeds = ne; +#endif + for(int i = 0; i < num_wave_speeds; ++i) { - max_char_speed = std::max(max_char_speed, ws[e]); + max_char_speed = std::max(max_char_speed, ws[i]); } return max_char_speed; @@ -178,8 +246,8 @@ namespace Theseus const int nfp = dc.num_face_points; const int nval_restr = operator_cache.restr_f->Height(); const int nfaces = nval_restr / (nfp * neq * 2); // (+/-) + const int npoints = nfaces * nfp; const int face_size = 2*nfp*neq; - const int norm_size = nfp*dim; if(operator_cache.uInt.Size() != nval_restr){ operator_cache.uInt.SetSize(nval_restr); @@ -211,11 +279,28 @@ namespace Theseus mfem::real_t *ws_d = dc.ifWaveSpeed_d; - mfem::forall(nfaces, [=] MFEM_HOST_DEVICE (int i) +#ifdef POINT_PARALLEL_INTERIOR_FACES + mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p) { - const int face_offset = i*face_size; - const int n_offset = i*norm_size; - const int w_offset = i*nfp; + const int f = p / nfp; + const int fp = p % nfp; + const int face_offset = f*face_size; + const int point_offset = f*nfp + fp; + + const mfem::real_t *u_face_d = u_d + face_offset; + mfem::real_t *rhs_face_d = rhs_d + face_offset; + + ws_d[p] = DGSEMIntegrator::AssembleFacePointKernel( + dc, u_face_d, nor_d + point_offset*dim, + inv1_d[point_offset], inv2_d[point_offset], fp, rhs_face_d); + }); +#else + const int norm_size = nfp*dim; + mfem::forall(nfaces, [=] MFEM_HOST_DEVICE (int f) + { + const int face_offset = f*face_size; + const int n_offset = f*norm_size; + const int w_offset = f*nfp; const mfem::real_t *u_face_d = u_d + face_offset; mfem::real_t *rhs_face_d = rhs_d + face_offset; @@ -225,9 +310,10 @@ namespace Theseus mfem::real_t ws = DGSEMIntegrator::AssembleElementFaceKernel(dc, u_face_d, nor_face_d, w_minus_d, w_plus_d, rhs_face_d); - ws_d[i] = ws; + ws_d[f] = ws; }); +#endif operator_cache.restr_f->MultTranspose(rhs_faces, faces_dudt); pdudt += faces_dudt; // on device? @@ -236,9 +322,14 @@ namespace Theseus // - Reduce for rank-local max_char_speed const mfem::real_t *ws = operator_cache.ifWaveSpeed.HostRead(); mfem::real_t max_char_speed_facial = 0.0; - for(int f = 0;f < operator_cache.num_interior_faces;f++) +#ifdef POINT_PARALLEL_INTERIOR_FACES + const int num_wave_speeds = npoints; +#else + const int num_wave_speeds = nfaces; +#endif + for(int i = 0; i < num_wave_speeds; ++i) { - max_char_speed_facial = std::max(max_char_speed_facial, ws[f]); + max_char_speed_facial = std::max(max_char_speed_facial, ws[i]); } return max_char_speed_facial; diff --git a/include/NSOperator_impl.hpp b/include/NSOperator_impl.hpp index d6772f5..9d9c165 100644 --- a/include/NSOperator_impl.hpp +++ b/include/NSOperator_impl.hpp @@ -237,6 +237,7 @@ namespace Theseus const int ne = dc.num_elements; const int ndof = dc.ndof_scalar_el; const int neq = dc.num_equations; + const int npoints = ne * ndof; const int estride = ndof * neq; const int jac_stride = ndof; const int metric_stride = ndof * dc.dim * dc.dim; @@ -244,6 +245,22 @@ namespace Theseus const mfem::real_t *elJac_d = dc.elJac_d; const mfem::real_t *elMetric_d = dc.elMetric_d; +#ifdef POINT_PARALLEL_VOLUME + mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p) + { + const int e = p / ndof; + const int point = p % ndof; + const int element_offset = e * estride; + mfem::real_t *du_el_d[Theseus::MAXDIM] = {nullptr, nullptr, nullptr}; + for(int idim = 0; idim < dim; ++idim){ + du_el_d[idim] = dU_d[idim] + element_offset; + } + + Theseus::DGSEMIntegrator::AssembleGradVolumePointKernel( + dc, Ue_d + element_offset, elJac_d + e*jac_stride, + elMetric_d + e*metric_stride, point, du_el_d); + }); +#else mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e) { const mfem::real_t *u_el = Ue_d + e * estride; @@ -258,6 +275,7 @@ namespace Theseus Theseus::DGSEMIntegrator::AssembleGradElementVolumeKernel(dc, u_el, jac_el, metric_el, du_el_d); }); +#endif for(int idim = 0;idim < dim;idim++){ operator_cache.restr_v->AddMultTranspose(dUe[idim], *p_grad_u[idim]); @@ -398,8 +416,8 @@ namespace Theseus const int neq = dc.num_equations; const int nfp = dc.num_face_points; const int nfaces = restr_size / (2 * nfp * neq); + const int npoints = nfaces * nfp; const int face_size = 2 * nfp * neq; - const int norm_size = nfp * dim; mfem::Vector &u_faces(operator_cache.sInt); if(u_faces.Size() != restr_size){ @@ -447,6 +465,25 @@ namespace Theseus const mfem::real_t *wm_d = dc.fw_minus_d; const mfem::real_t *wp_d = dc.fw_plus_d; +#ifdef POINT_PARALLEL_INTERIOR_FACES + mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p) + { + const int f = p / nfp; + const int fp = p % nfp; + const int face_offset = f * face_size; + const int point_offset = f * nfp + fp; + + mfem::real_t *rhs_face[Theseus::MAXDIM] = {nullptr, nullptr, nullptr}; + for(int idim = 0; idim < dim; ++idim){ + rhs_face[idim] = rhs_d[idim] + face_offset; + } + + Theseus::DGSEMIntegrator::AssembleGradInteriorFacePointKernel( + dc, u_d + face_offset, nor_d + point_offset*dim, + wm_d[point_offset], wp_d[point_offset], fp, rhs_face); + }); +#else + const int norm_size = nfp * dim; mfem::forall(nfaces, [=] MFEM_HOST_DEVICE (int f) { const int face_offset = f * face_size; @@ -470,6 +507,7 @@ namespace Theseus w_plus_d, rhs_face); }); +#endif for(int idim = 0;idim < dim;idim++){ operator_cache.restr_f->MultTranspose(rhs_faces[idim], duInt); @@ -550,8 +588,8 @@ namespace Theseus const int neq = dc.num_equations; const int nfp = dc.num_face_points; const int nfaces = operator_cache.restr_f->Height() / (nfp * neq * 2); // (+/-) + const int npoints = nfaces * nfp; const int face_size = 2*nfp*neq; - const int norm_size = nfp*dim; const int restr_size = operator_cache.restr_f->Height(); mfem::Vector &int_u(operator_cache.uInt); @@ -610,11 +648,37 @@ namespace Theseus mfem::real_t *ws_d = dc.ifWaveSpeed_d; - mfem::forall(nfaces, [=] MFEM_HOST_DEVICE (int i) +#ifdef POINT_PARALLEL_INTERIOR_FACES + mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p) + { + const int f = p / nfp; + const int fp = p % nfp; + const int face_offset = f*face_size; + const int point_offset = f*nfp + fp; + + const mfem::real_t *u_face_d = u_d + face_offset; + mfem::real_t *rhs_face_d = rhs_d + face_offset; + const mfem::real_t *dprim_face_x = + (dim > 0) ? grad_prim_d[0] + face_offset : nullptr; + const mfem::real_t *dprim_face_y = + (dim > 1) ? grad_prim_d[1] + face_offset : nullptr; + const mfem::real_t *dprim_face_z = + (dim > 2) ? grad_prim_d[2] + face_offset : nullptr; + const mfem::real_t radius = + dc.axisymmetric ? face_radius_d[point_offset] : 0.0; + + ws_d[p] = Theseus::DGSEMIntegrator::AssembleViscousFacePointKernel( + dc, u_face_d, nor_d + point_offset*dim, + inv1_d[point_offset], inv2_d[point_offset], + dprim_face_x, dprim_face_y, dprim_face_z, radius, fp, rhs_face_d); + }); +#else + const int norm_size = nfp*dim; + mfem::forall(nfaces, [=] MFEM_HOST_DEVICE (int f) { - const int face_offset = i*face_size; - const int n_offset = i*norm_size; - const int w_offset = i*nfp; + const int face_offset = f*face_size; + const int n_offset = f*norm_size; + const int w_offset = f*nfp; const mfem::real_t *u_face_d = u_d + face_offset; mfem::real_t *rhs_face_d = rhs_d + face_offset; @@ -635,8 +699,9 @@ namespace Theseus dprim_face_z, radius_face_d, rhs_face_d); - ws_d[i] = ws; + ws_d[f] = ws; }); +#endif operator_cache.restr_f->MultTranspose(rhs_faces, faces_dudt); pdudt += faces_dudt; // on device? @@ -645,9 +710,14 @@ namespace Theseus // - Reduce for rank-local max_char_speed const mfem::real_t *ws = operator_cache.ifWaveSpeed.HostRead(); mfem::real_t max_char_speed_facial = 0.0; - for(int f = 0;f < operator_cache.num_interior_faces;f++) +#ifdef POINT_PARALLEL_INTERIOR_FACES + const int num_wave_speeds = npoints; +#else + const int num_wave_speeds = nfaces; +#endif + for(int i = 0; i < num_wave_speeds; ++i) { - max_char_speed_facial = std::max(max_char_speed_facial, ws[f]); + max_char_speed_facial = std::max(max_char_speed_facial, ws[i]); } return max_char_speed_facial; @@ -895,6 +965,25 @@ namespace Theseus mfem::real_t *ws_d = dc.elWaveSpeed_d; +#ifdef POINT_PARALLEL_VOLUME + const int npoints = ne * ndof; + mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p) + { + const int e = p / ndof; + const int point = p % ndof; + const int attr = elem_attr_d[e]; + if (attr_marker_d[attr-1] == 0) { + ws_d[p] = 0.0; + return; + } + + const int element_offset = e * estride; + ws_d[p] = DGSEMIntegrator::AssembleVolumePointKernel( + dc, Ue_d + element_offset, elJac_d + e*jac_stride, + elMetric_d + e*metric_stride, point, dUe_d + element_offset); + }); +#endif + // Inside the FORALL below, executed on device mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e) { @@ -915,9 +1004,13 @@ namespace Theseus const mfem::real_t *u_el = Ue_d + eoff; mfem::real_t *du_el = dUe_d + eoff; - mfem::real_t cs_el = \ - Theseus::DGSEMIntegrator::AssembleElementVolumeKernel(dc, u_el, - jac_el, metric_el, du_el); +#ifdef POINT_PARALLEL_VOLUME + mfem::real_t cs_el = ws_d[e*ndof]; +#else + mfem::real_t cs_el = + Theseus::DGSEMIntegrator::AssembleElementVolumeKernel( + dc, u_el, jac_el, metric_el, du_el); +#endif #ifdef SUBCELL_FV_BLENDING mfem::real_t alpha_fv = alpha_d[e]; if(alpha_fv > 1e-16){ @@ -940,11 +1033,13 @@ namespace Theseus #endif AddAxisymmetricEulerElementSource( dc, u_el, radius_el, jac_el, metric_el, du_el); +#ifdef POINT_PARALLEL_VOLUME + ws_d[e*ndof] = cs_el; +#else ws_d[e] = cs_el; +#endif - // Inviscid part is done: dUe currrently holds the inviscid RHS - // Host code mixes inviscid and viscous assembly, we need separate. - // Call the Viscous Assembly routine +#ifndef POINT_PARALLEL_VOLUME const mfem::real_t *grad_prim_el[Theseus::MAXDIM] = {nullptr, nullptr, nullptr}; for(int idim = 0;idim < dim;idim++){ grad_prim_el[idim] = gradPrim_d[idim] + eoff; @@ -954,8 +1049,35 @@ namespace Theseus radius_el, grad_prim_el[0], grad_prim_el[1], grad_prim_el[2], du_el); +#endif }); + +#ifdef POINT_PARALLEL_VOLUME + mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p) + { + const int e = p / ndof; + const int point = p % ndof; + const int attr = elem_attr_d[e]; + if (attr_marker_d[attr-1] == 0) { + return; + } + + const int element_offset = e * estride; + const mfem::real_t *grad_prim_el[Theseus::MAXDIM] = { + nullptr, nullptr, nullptr}; + for(int idim = 0; idim < dim; ++idim){ + grad_prim_el[idim] = gradPrim_d[idim] + element_offset; + } + + Theseus::DGSEMIntegrator::AssembleViscousVolumePointKernel( + dc, Ue_d + element_offset, elJac_d + e*jac_stride, + elMetric_d + e*metric_stride, + dc.axisymmetric ? elRadius_d + e*jac_stride : nullptr, + grad_prim_el[0], grad_prim_el[1], grad_prim_el[2], point, + dUe_d + element_offset); + }); +#endif // The rest is identical to Euler operator // Scatter RHS back to storage @@ -965,9 +1087,14 @@ namespace Theseus // - Reduce for rank-local max_char_speed const mfem::real_t *ws = operator_cache.elWaveSpeed.HostRead(); mfem::real_t max_char_speed = 0.0; - for(int e = 0;e < operator_cache.num_elements;e++) +#ifdef POINT_PARALLEL_VOLUME + const int num_wave_speeds = ne * ndof; +#else + const int num_wave_speeds = ne; +#endif + for(int i = 0; i < num_wave_speeds; ++i) { - max_char_speed = std::max(max_char_speed, ws[e]); + max_char_speed = std::max(max_char_speed, ws[i]); } return max_char_speed; diff --git a/include/dgsem_cache.hpp b/include/dgsem_cache.hpp index a87b816..88ca199 100644 --- a/include/dgsem_cache.hpp +++ b/include/dgsem_cache.hpp @@ -120,8 +120,9 @@ namespace Theseus mfem::Vector bc_vector_data; // Physics parts - used directly on device - mutable mfem::Vector elWaveSpeed; // size nelements - mutable mfem::Vector ifWaveSpeed; // size ninterior faces + // One value per element, or per element point in point-volume builds. + mutable mfem::Vector elWaveSpeed; + mutable mfem::Vector ifWaveSpeed; // size ninterior faces * points per face mutable mfem::Vector bndWaveSpeed; // size nbnd faces OperatorGasModel gas; InviscidFlux iflux; diff --git a/include/dgsem_cache_utilities.hpp b/include/dgsem_cache_utilities.hpp index 78e1648..44e511e 100644 --- a/include/dgsem_cache_utilities.hpp +++ b/include/dgsem_cache_utilities.hpp @@ -171,7 +171,11 @@ namespace Theseus { std::memcpy(cache->Dhat.HostWrite(), Dhat_T.Data(), sizeof(mfem::real_t)*Np_x*Np_x); std::memcpy(cache->Dhat2.HostWrite(), Dhat2_T.Data(), sizeof(mfem::real_t)*Np_x*Np_x); +#ifdef POINT_PARALLEL_VOLUME + cache->elWaveSpeed.SetSize(nelem * cache->ndof_scalar_el); +#else cache->elWaveSpeed.SetSize(nelem); +#endif cache->elWaveSpeed = 0.0; cache->elWaveSpeed.UseDevice(); cache->elWaveSpeed.Read(); @@ -208,7 +212,7 @@ namespace Theseus { cache->face_wt_plus.Read(); cache->face_radius.Read(); - cache->ifWaveSpeed.SetSize(cache->num_interior_faces); + cache->ifWaveSpeed.SetSize(cache->num_interior_faces * nfp); cache->ifWaveSpeed = 0.0; cache->ifWaveSpeed.UseDevice(); cache->ifWaveSpeed.Read(); @@ -1045,7 +1049,15 @@ namespace Theseus { MFEM_VERIFY(ds_size > 0, "Elem attr not set"); ds_size = cache.elWaveSpeed.Size(); +#ifdef POINT_PARALLEL_VOLUME + MFEM_VERIFY(ds_size == cache.num_elements * cache.ndof_scalar_el, + "Element point wavespeeds missized."); +#else MFEM_VERIFY(ds_size == cache.num_elements, "Element wavespeeds missized."); +#endif + ds_size = cache.ifWaveSpeed.Size(); + MFEM_VERIFY(ds_size == cache.num_interior_faces * cache.num_face_points, + "Interior-face wavespeeds missized."); ds_size = cache.bndWaveSpeed.Size(); ds_size = cache.elJac.Size(); MFEM_VERIFY(ds_size > 0, "Element Jacobians not set"); diff --git a/include/timer.hpp b/include/timer.hpp index c9c4fa1..dfe894b 100644 --- a/include/timer.hpp +++ b/include/timer.hpp @@ -10,6 +10,11 @@ #ifdef TIMER_SYNC_DEVICE #include "mfem.hpp" #endif +#include +#include +#include +#include +#include namespace Theseus { @@ -27,7 +32,7 @@ namespace Theseus : name_(name) { MFEM_DEVICE_SYNC; - // mfem::Device::Sync(); + // mfem::Device::Sync(); start_ = clock::now(); } #endif @@ -49,7 +54,7 @@ namespace Theseus MPI_Barrier(MPI_COMM_WORLD); end = clock::now(); if(rank == 0) - MPI_Reduce(MPI_IN_PLACE,&global_ms, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + MPI_Reduce(MPI_IN_PLACE,&global_ms, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); else MPI_Reduce(&global_ms, NULL, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); #endif @@ -77,4 +82,258 @@ namespace Theseus clock::time_point start_; }; + + class TimestepTimer + { + public: + using clock = std::chrono::steady_clock; + + explicit TimestepTimer(MPI_Comm comm = MPI_COMM_WORLD) + : comm_(comm) + { + MPI_Comm_rank(comm_, &rank_); + MPI_Comm_size(comm_, &nranks_); + } + + void Start() + { + Sync(); + start_ = clock::now(); + running_ = true; + } + + void Stop() + { + if (!running_) + { + return; + } + + Sync(); + const auto stop = clock::now(); + + const double elapsed_ms = + std::chrono::duration(stop - start_).count(); + + total_ms_ += elapsed_ms; + min_ms_ = std::min(min_ms_, elapsed_ms); + max_ms_ = std::max(max_ms_, elapsed_ms); + ++count_; + + running_ = false; + } + + void Reset() + { + count_ = 0; + total_ms_ = 0.0; + min_ms_ = std::numeric_limits::max(); + max_ms_ = 0.0; + running_ = false; + } + + std::uint64_t Count() const + { + return count_; + } + + double Total() const + { + return total_ms_; + } + + double Min() const + { + return count_ > 0 ? min_ms_ : 0.0; + } + + double Max() const + { + return count_ > 0 ? max_ms_ : 0.0; + } + + double Mean() const + { + return count_ > 0 ? total_ms_ / static_cast(count_) : 0.0; + } + + void Finalize(std::ostream &os = std::cout) const + { + // + // Local statistics. + // + const double local_total = Total(); + const double local_min = Min(); + const double local_max = Max(); + const double local_mean = Mean(); + + const unsigned long long local_count = + static_cast(count_); + + // + // Count statistics across ranks. + // + unsigned long long count_min = 0; + unsigned long long count_max = 0; + unsigned long long count_sum = 0; + + MPI_Reduce(&local_count, &count_min, 1, + MPI_UNSIGNED_LONG_LONG, MPI_MIN, 0, comm_); + + MPI_Reduce(&local_count, &count_max, 1, + MPI_UNSIGNED_LONG_LONG, MPI_MAX, 0, comm_); + + MPI_Reduce(&local_count, &count_sum, 1, + MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, comm_); + + // + // For each local statistic, collect min/max/sum across ranks. + // + double total_min = 0.0; + double total_max = 0.0; + double total_sum = 0.0; + + double step_min_min = 0.0; + double step_min_max = 0.0; + double step_min_sum = 0.0; + + double step_max_min = 0.0; + double step_max_max = 0.0; + double step_max_sum = 0.0; + + double mean_min = 0.0; + double mean_max = 0.0; + double mean_sum = 0.0; + + ReduceMinMaxSum(local_total, + total_min, total_max, total_sum); + + ReduceMinMaxSum(local_min, + step_min_min, step_min_max, step_min_sum); + + ReduceMinMaxSum(local_max, + step_max_min, step_max_max, step_max_sum); + + ReduceMinMaxSum(local_mean, + mean_min, mean_max, mean_sum); + + if (rank_ != 0) + { + return; + } + + const double inv_nranks = 1.0 / static_cast(nranks_); + + const double count_mean = + static_cast(count_sum) * inv_nranks; + + const double total_mean = + total_sum * inv_nranks; + + const double step_min_mean = + step_min_sum * inv_nranks; + + const double step_max_mean = + step_max_sum * inv_nranks; + + const double mean_mean = + mean_sum * inv_nranks; + + os << std::fixed << std::setprecision(3); + + os << "\n" + << "========================================================================\n" + << "Theseus timestep performance\n" + << "========================================================================\n" + << "MPI ranks : " << nranks_ << "\n" + << "\n" + << "Per-rank timestep statistics\n" + << " min" + << " mean" + << " max\n" + << " Timed steps : " + << std::setw(12) << count_min + << std::setw(14) << count_mean + << std::setw(14) << count_max << "\n" + << " Total time (ms) : " + << std::setw(12) << total_min + << std::setw(14) << total_mean + << std::setw(14) << total_max << "\n" + << " Min step (ms) : " + << std::setw(12) << step_min_min + << std::setw(14) << step_min_mean + << std::setw(14) << step_min_max << "\n" + << " Mean step (ms) : " + << std::setw(12) << mean_min + << std::setw(14) << mean_mean + << std::setw(14) << mean_max << "\n" + << " Max step (ms) : " + << std::setw(12) << step_max_min + << std::setw(14) << step_max_mean + << std::setw(14) << step_max_max << "\n" + << "\n" + << "Critical-rank performance\n" + << " Mean timestep : " << mean_max << " ms\n"; + + if (mean_max > 0.0) + { + os << " Timesteps/sec : " + << 1000.0 / mean_max << "\n"; + } + +#ifdef TIMER_SYNC_DEVICE + os << " Device sync : enabled\n"; +#else + os << " Device sync : disabled\n"; +#endif + +#ifdef TIMER_BARRIER + os << " MPI barrier : enabled\n"; +#else + os << " MPI barrier : disabled\n"; +#endif + + os << "========================================================================\n"; + } + + private: + void Sync() const + { +#ifdef TIMER_BARRIER + MPI_Barrier(comm_); +#endif + +#ifdef TIMER_SYNC_DEVICE + MFEM_DEVICE_SYNC; +#endif + } + + void ReduceMinMaxSum(double local, + double &minimum, + double &maximum, + double &sum) const + { + MPI_Reduce(&local, &minimum, 1, + MPI_DOUBLE, MPI_MIN, 0, comm_); + + MPI_Reduce(&local, &maximum, 1, + MPI_DOUBLE, MPI_MAX, 0, comm_); + + MPI_Reduce(&local, &sum, 1, + MPI_DOUBLE, MPI_SUM, 0, comm_); + } + + MPI_Comm comm_ = MPI_COMM_WORLD; + int rank_ = 0; + int nranks_ = 1; + + clock::time_point start_; + + std::uint64_t count_ = 0; + double total_ms_ = 0.0; + double min_ms_ = std::numeric_limits::max(); + double max_ms_ = 0.0; + + bool running_ = false; + }; } diff --git a/src/Simulation.cpp b/src/Simulation.cpp index bb06f7a..351f5ef 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -322,17 +322,17 @@ namespace Theseus if (exact_signature == 0) { exact_solution = std::make_unique( - num_equations, - Prandtl::ConditionFactory::Instance(). - GetVectorTDFunctionBoundaryCondition0(key)()); + num_equations, + Prandtl::ConditionFactory::Instance(). + GetVectorTDFunctionBoundaryCondition0(key)()); } else if (exact_signature == 1) { const mfem::real_t x1 = exact["params"].value("x1", 0.0); exact_solution = std::make_unique( - num_equations, - Prandtl::ConditionFactory::Instance(). - GetVectorTDFunctionBoundaryCondition1(key)(x1)); + num_equations, + Prandtl::ConditionFactory::Instance(). + GetVectorTDFunctionBoundaryCondition1(key)(x1)); } else { @@ -471,13 +471,13 @@ namespace Theseus std::int64_t num_elements_local = ndofscalar / points_per_element; std::int64_t num_elements_total = num_elements_local; -if (debug_simulation && numProcs > 1) { + if (debug_simulation && numProcs > 1) { for(int irank = 0;irank < numProcs;irank++){ - if(myRank == irank){ - std::cout << "Rank(" << myRank << ") Number of elements: " - << num_elements_local << std::endl; - } - MPI_Barrier(pmesh->GetComm()); + if(myRank == irank){ + std::cout << "Rank(" << myRank << ") Number of elements: " + << num_elements_local << std::endl; + } + MPI_Barrier(pmesh->GetComm()); } } if(numProcs > 1){ @@ -486,8 +486,8 @@ if (debug_simulation && numProcs > 1) { MPI_Allreduce(MPI_IN_PLACE, &max_nel, 1, MPI_LONG_LONG, MPI_MAX, pmesh->GetComm()); MPI_Allreduce(MPI_IN_PLACE, &min_nel, 1, MPI_LONG_LONG, MPI_MIN, pmesh->GetComm()); if(myRank == 0){ - std::cout << "Partition NumElements (min, max) = (" << min_nel << ", " << max_nel - << ")" << std::endl; + std::cout << "Partition NumElements (min, max) = (" << min_nel << ", " << max_nel + << ")" << std::endl; } } MPI_Allreduce(MPI_IN_PLACE, &num_elements_total, 1, MPI_LONG_LONG, MPI_SUM, pmesh->GetComm()); @@ -754,38 +754,38 @@ if (debug_simulation && numProcs > 1) { { bc_descr.type = int(Theseus::BCType::NoSlipIso); bc_descr.data_kind = int(Theseus::BCDataKind::VectorAndScalarConstant); -if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") && - bc_props.contains("temperature") && bc_props["temperature"].contains("scalar"))) - { std::cerr << "Error: no-slip-isothermal requires velocity.vector and temperature.scalar." << std::endl; return 1; } -{ - std::string velBC_key = bc_props["velocity"]["vector"].get(); - std::string tempBC_key = bc_props["temperature"]["scalar"].get(); - // std::string state_key = bc_props["vector"].get(); - auto vel_bc = Prandtl::ConditionFactory::Instance().GetVectorBoundaryCondition(velBC_key); - auto temp_bc = Prandtl::ConditionFactory::Instance().GetScalarBoundaryCondition(tempBC_key); - - mfem::Vector bc_data(vel_bc.Size() + 1); - std::ostringstream Ostr; - Ostr << "Wall velocity: < "; - for(int ivec = 0;ivec < vel_bc.Size();ivec++){ - bc_data[ivec] = vel_bc[ivec]; - Ostr << vel_bc[ivec] << " "; - } - Ostr << ">" << std::endl; - Ostr << "Wall temperature: " << temp_bc << std::endl; - bc_data[vel_bc.Size()] = temp_bc; - bc_descr.data_index = Theseus::AppendBCVectorPayload(bc_vector_data, bc_data); - Ostr << "bc_vector_data index: " << bc_descr.data_index << std::endl - << "BC Data So Far: ["; - for(int ivec=0;ivec < bc_vector_data.Size();ivec++){ - Ostr << bc_vector_data[ivec] << " "; - } - Ostr << "]" << std::endl; - if(debug_simulation && mfem::Mpi::Root()){ - std::cout << Ostr.str(); - } - rhsOp->AddBdrFaceMarker(bdr_marker_vector.back()); - } + if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") && + bc_props.contains("temperature") && bc_props["temperature"].contains("scalar"))) + { std::cerr << "Error: no-slip-isothermal requires velocity.vector and temperature.scalar." << std::endl; return 1; } + { + std::string velBC_key = bc_props["velocity"]["vector"].get(); + std::string tempBC_key = bc_props["temperature"]["scalar"].get(); + // std::string state_key = bc_props["vector"].get(); + auto vel_bc = Prandtl::ConditionFactory::Instance().GetVectorBoundaryCondition(velBC_key); + auto temp_bc = Prandtl::ConditionFactory::Instance().GetScalarBoundaryCondition(tempBC_key); + + mfem::Vector bc_data(vel_bc.Size() + 1); + std::ostringstream Ostr; + Ostr << "Wall velocity: < "; + for(int ivec = 0;ivec < vel_bc.Size();ivec++){ + bc_data[ivec] = vel_bc[ivec]; + Ostr << vel_bc[ivec] << " "; + } + Ostr << ">" << std::endl; + Ostr << "Wall temperature: " << temp_bc << std::endl; + bc_data[vel_bc.Size()] = temp_bc; + bc_descr.data_index = Theseus::AppendBCVectorPayload(bc_vector_data, bc_data); + Ostr << "bc_vector_data index: " << bc_descr.data_index << std::endl + << "BC Data So Far: ["; + for(int ivec=0;ivec < bc_vector_data.Size();ivec++){ + Ostr << bc_vector_data[ivec] << " "; + } + Ostr << "]" << std::endl; + if(debug_simulation && mfem::Mpi::Root()){ + std::cout << Ostr.str(); + } + rhsOp->AddBdrFaceMarker(bdr_marker_vector.back()); + } } else if (type == "supersonic-outflow") { @@ -1043,7 +1043,7 @@ if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") & pd->SetLevelsOfDetail(order); pd->SetDataFormat(mfem::VTKFormat::BINARY); pd->SetHighOrderOutput( - visualization_config.MeshMode() == VisualizationMeshMode::vtk_high_order); + visualization_config.MeshMode() == VisualizationMeshMode::vtk_high_order); } else if (visit) { @@ -1084,7 +1084,7 @@ if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") & if (mfem::Mpi::Root()) { std::cout << "================================================" << std::endl - << "Theseus Simulation Running Now" << std::endl + << "Theseus::Simulation Running " << std::endl << "================================================" << std::endl; } const auto &gasModel = rhsOp->GetGasModelInterface(); @@ -1206,154 +1206,175 @@ if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") & if(mfem::Mpi::Root()){ std::cout << "Writing initial soln..." << std::endl; } - Theseus::ScopedTimer timer("VisInit"); + { + Theseus::ScopedTimer timer("VisInit"); - UpdateVisualizationFields(); + UpdateVisualizationFields(); - SaveVisualization(); + SaveVisualization(); + } + if(mfem::Mpi::Root()){ + std::cout << "Done writing initial soln." << std::endl; + } } - while (!done) - { - - if (debug_simulation) - { - MPI_Barrier(pmesh->GetComm()); - if(mfem::Mpi::Root()){ - std::cout << "############################################" - << std::endl - << "[TIME STEP = " << ti << ", TIME = " << t << "]" - << std::endl - << "############################################" - << std::endl; - } - } - - // Compute the time step size - dt_real = std::min(dt, t_final - t); - - // Perform the time step - { - Theseus::ScopedTimer timer("Timestep"); - ode_solver->Step(*sol, t, dt_real); - } - ti++; - - mfem::real_t cfl_rep = 0.0; - if (ti % print_interval == 0 || (variable_dt && cfl > 0) || debug_simulation){ - rhsOp->ComputeIntegralMeasures(*sol, diag); - } - // Update the time step size with CFL? - if ((variable_dt && cfl > 0) || (ti%print_interval == 0) || debug_simulation) - { - mfem::real_t max_char_speed = rhsOp->GetMaxCharSpeed(); - MPI_Allreduce(MPI_IN_PLACE, &max_char_speed, 1, mfem::MPITypeMap::mpi_type, - MPI_MAX, pmesh->GetComm()); - mfem::real_t dt_adv = heff / max_char_speed; - mfem::real_t dt_est = dt_adv; + if(mfem::Mpi::Root()){ + std::cout << "Timestepping ..." << std::endl; + } + { + TimestepTimer timestep_timer; + ScopedTimer timestepping_timer("Timestepping"); + while (!done) + { + + if (debug_simulation) + { + MPI_Barrier(pmesh->GetComm()); + if(mfem::Mpi::Root()){ + std::cout << "############################################" + << std::endl + << "[TIME STEP = " << ti << ", TIME = " << t << "]" + << std::endl + << "############################################" + << std::endl; + } + } + + // Compute the time step size + dt_real = std::min(dt, t_final - t); + + // Perform the time step + if(ti > 1){ + timestep_timer.Start(); + } + { + Theseus::ScopedTimer timer("Timestep"); + ode_solver->Step(*sol, t, dt_real); + } + if(ti > 1){ + timestep_timer.Stop(); + } + ti++; + + mfem::real_t cfl_rep = 0.0; + if (ti % print_interval == 0 || (variable_dt && cfl > 0) || debug_simulation){ + rhsOp->ComputeIntegralMeasures(*sol, diag); + } + // Update the time step size with CFL? + if ((variable_dt && cfl > 0) || (ti%print_interval == 0) || debug_simulation) + { + mfem::real_t max_char_speed = rhsOp->GetMaxCharSpeed(); + MPI_Allreduce(MPI_IN_PLACE, &max_char_speed, 1, mfem::MPITypeMap::mpi_type, + MPI_MAX, pmesh->GetComm()); + mfem::real_t dt_adv = heff / max_char_speed; + mfem::real_t dt_est = dt_adv; #ifdef PARABOLIC - mfem::real_t nu_eff = nuscale * physicsConstants.mu / diag.min_dens; - mfem::real_t dt_diff = heff * heff / nu_eff; - mfem::real_t dt_m1 = 1.0 / (1.0/dt_adv + 1.0/dt_diff); - dt_est = dt_m1; + mfem::real_t nu_eff = nuscale * physicsConstants.mu / diag.min_dens; + mfem::real_t dt_diff = heff * heff / nu_eff; + mfem::real_t dt_m1 = 1.0 / (1.0/dt_adv + 1.0/dt_diff); + dt_est = dt_m1; #endif - if(variable_dt){ - dt = cfl / dim * dt_est; - } else { - cfl_rep = dim * dt / dt_est; - } + if(variable_dt){ + dt = cfl / dim * dt_est; + } else { + cfl_rep = dim * dt / dt_est; + } - if(debug_simulation && mfem::Mpi::Root()){ + if(debug_simulation && mfem::Mpi::Root()){ #ifdef PARABOLIC - std::cout << "DT(adv, diff, sim): (" << dt_adv << ", " << dt_diff - << ", " << dt << ")" << std::endl - << "Effective viscosity: " << nu_eff << std::endl; + std::cout << "DT(adv, diff, sim): (" << dt_adv << ", " << dt_diff + << ", " << dt << ")" << std::endl + << "Effective viscosity: " << nu_eff << std::endl; #else - std::cout << "DT(adv, sim): (" << dt_adv << ", " << dt << ")" << std::endl; + std::cout << "DT(adv, sim): (" << dt_adv << ", " << dt << ")" << std::endl; #endif - if(!variable_dt){ - std::cout << "CFL: "<< cfl_rep << std::endl; - } - std::cout << "Max wavespeed: " << max_char_speed << std::endl - << "Max specific volume: " << 1.0 / diag.min_dens << std::endl; - } - - } - - // Check for completion - done = ((t >= t_final - 1e-8 * dt) || StepLimitReached(ti, nsteps_max)); - - // Check for NaN/Inf values? - rho.HostRead(); - if (nancheck && ti % nancheck_steps == 0) - { - for (const mfem::real_t &val : rho) - { - if (std::isnan(val) || std::isinf(val)) - { - MFEM_ABORT("NaN/Inf Detected at Time Step " + std::to_string(ti) + - " on Rank " + std::to_string(myRank)); - break; - } - } - } - // Visualize the solution? - if (visualize && (done || t >= next_save_t || ti % vis_steps == 0)) - { - - UpdateVisualizationFields(); - - SaveVisualization(); - - - save_dt = (t < trigger_t) ? save_dt1 : save_dt2; - next_save_t += save_dt; - - } - - - if (checkpoint_config.SaveEnabled() && (done || t >= next_checkpoint_t)) - { - SaveCheckpoint(); - next_checkpoint_t += checkpoint_config.Interval(); - } - - if (ti % print_interval == 0 || debug_simulation) - { - mfem::real_t ke0 = diag0.ke; - if(ke0 == 0.0){ ke0 = 1.0; }; - if (mfem::Mpi::Root()) - { - std::ostringstream Ostr; - Ostr << "time step: " << ti << ", time: " << t; - if(variable_dt){ - Ostr << ", dt: " << dt; - } else { - Ostr << ", cfl: " << cfl_rep; - } - Ostr << std::endl - << "rho(" << diag.min_dens << "," << diag.max_dens << "), " - << "p(" << diag.min_press << "," << diag.max_press << "), " - << "T(" << diag.min_temp << "," << diag.max_temp << ")" << std::endl - << "TotalChange: Mass: " << (diag.mass - diag0.mass) / diag0.mass - << ", Energy: " << (diag.en - diag0.en) / diag0.en - << ", K.E.: " << (diag.ke - diag0.ke) / ke0 << std::endl; - std::cout << Ostr.str(); - } - } - } - + if(!variable_dt){ + std::cout << "CFL: "<< cfl_rep << std::endl; + } + std::cout << "Max wavespeed: " << max_char_speed << std::endl + << "Max specific volume: " << 1.0 / diag.min_dens << std::endl; + } + + } + + // Check for completion + done = ((t >= t_final - 1e-8 * dt) || StepLimitReached(ti, nsteps_max)); + + // Check for NaN/Inf values? + rho.HostRead(); + if (nancheck && ti % nancheck_steps == 0) + { + for (const mfem::real_t &val : rho) + { + if (std::isnan(val) || std::isinf(val)) + { + MFEM_ABORT("NaN/Inf Detected at Time Step " + std::to_string(ti) + + " on Rank " + std::to_string(myRank)); + break; + } + } + } + // Visualize the solution? + if (visualize && (done || t >= next_save_t || ti % vis_steps == 0)) + { + + UpdateVisualizationFields(); + + SaveVisualization(); + + + save_dt = (t < trigger_t) ? save_dt1 : save_dt2; + next_save_t += save_dt; + + } + + + if (checkpoint_config.SaveEnabled() && (done || t >= next_checkpoint_t)) + { + SaveCheckpoint(); + next_checkpoint_t += checkpoint_config.Interval(); + } + + if (ti % print_interval == 0 || debug_simulation) + { + mfem::real_t ke0 = diag0.ke; + if(ke0 == 0.0){ ke0 = 1.0; }; + if (mfem::Mpi::Root()) + { + std::ostringstream Ostr; + Ostr << "time step: " << ti << ", time: " << t; + if(variable_dt){ + Ostr << ", dt: " << dt; + } else { + Ostr << ", cfl: " << cfl_rep; + } + Ostr << std::endl + << "rho(" << diag.min_dens << "," << diag.max_dens << "), " + << "p(" << diag.min_press << "," << diag.max_press << "), " + << "T(" << diag.min_temp << "," << diag.max_temp << ")" << std::endl + << "TotalChange: Mass: " << (diag.mass - diag0.mass) / diag0.mass + << ", Energy: " << (diag.en - diag0.en) / diag0.en + << ", K.E.: " << (diag.ke - diag0.ke) / ke0 << std::endl; + std::cout << Ostr.str(); + } + } + } + timestep_timer.Finalize(); + } + if(mfem::Mpi::Root()){ + std::cout << "Timestepping done." << std::endl; + } if (exact_solution) { exact_solution->SetTime(t); mfem::FunctionCoefficient cylindrical_weight( - [](const mfem::Vector &x) { - return AxisymmetricGeometry::MeasureMultiplier( - AxisymmetricGeometry::enabled, - AxisymmetricGeometry::enabled ? - AxisymmetricGeometry::Radius(x.GetData()) : 0.0); - }); + [](const mfem::Vector &x) { + return AxisymmetricGeometry::MeasureMultiplier( + AxisymmetricGeometry::enabled, + AxisymmetricGeometry::enabled ? + AxisymmetricGeometry::Radius(x.GetData()) : 0.0); + }); const mfem::real_t l1_error = sol->ComputeLpError(1.0, *exact_solution, &cylindrical_weight); const mfem::real_t l2_error = @@ -1458,13 +1479,13 @@ if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") & constexpr bool axisymmetric = false; #endif return {numProcs, - order, - dim, - num_equations, - static_cast(sizeof(mfem::real_t)), - static_cast(pmesh->GetGlobalNE()), - static_cast(vfes->GlobalVSize()), - axisymmetric}; + order, + dim, + num_equations, + static_cast(sizeof(mfem::real_t)), + static_cast(pmesh->GetGlobalNE()), + static_cast(vfes->GlobalVSize()), + axisymmetric}; } void Simulation::LoadCheckpoint() @@ -1487,7 +1508,7 @@ if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") & try { current_format = CheckpointConfig::ValidateMetadata( - metadata, CurrentCheckpointCompatibility()); + metadata, CurrentCheckpointCompatibility()); if (!current_format && mfem::Mpi::Root()) { std::cerr << "Warning: loading a legacy checkpoint without compatibility metadata; " @@ -1520,8 +1541,8 @@ if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") & << " does not match current finite-element space size " << vfes->GetVSize()); sol = std::make_shared(vfes.get()); checkpoint_stream.read( - reinterpret_cast(sol->HostWrite()), - static_cast(stored_size * sizeof(mfem::real_t))); + reinterpret_cast(sol->HostWrite()), + static_cast(stored_size * sizeof(mfem::real_t))); MFEM_VERIFY(checkpoint_stream, "Failed while reading checkpoint state: " << rank_file); } @@ -1569,8 +1590,8 @@ if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") & checkpoint_stream << "THESEUS_CHECKPOINT_RAW_V1\n"; checkpoint_stream.write(reinterpret_cast(&state_size), sizeof(state_size)); checkpoint_stream.write( - reinterpret_cast(sol->HostRead()), - static_cast(state_size * sizeof(mfem::real_t))); + reinterpret_cast(sol->HostRead()), + static_cast(state_size * sizeof(mfem::real_t))); checkpoint_stream.close(); MFEM_VERIFY(checkpoint_stream, "Failed while writing checkpoint state: " << rank_file); @@ -1594,6 +1615,4 @@ if (!(bc_props.contains("velocity") && bc_props["velocity"].contains("vector") & } MPI_Barrier(pmesh->GetComm()); } - - }