diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index 89684e997b..13d52326c3 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -574,6 +574,17 @@ public: /// amrex::Real estTimeStep (amrex::Real dt_old); + +/// +/// Compute the CFL timestep +/// + amrex::Real estdt_cfl(const amrex::Real time); + +/// +/// Diffusion-limited timestep +/// + amrex::Real estdt_temp_diffusion(); + /// /// Compute initial time step. /// diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 891a0f0614..5583132dd2 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -1287,7 +1287,7 @@ Castro::estTimeStep (Real dt_old) Real estdt = max_dt; const MultiFab& stateMF = get_new_data(State_Type); - + Real time = state[State_Type].curTime(); const Real* dx = geom.CellSize(); std::string limiter = "castro.max_dt"; @@ -1337,24 +1337,7 @@ Castro::estTimeStep (Real dt_old) { #endif -#ifdef _OPENMP -#pragma omp parallel reduction(min:estdt_hydro) -#endif - { - Real dt = max_dt / cfl; - - for (MFIter mfi(stateMF, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& box = mfi.tilebox(); - -#pragma gpu box(box) - ca_estdt(AMREX_INT_ANYD(box.loVect()), AMREX_INT_ANYD(box.hiVect()), - BL_TO_FORTRAN_ANYD(stateMF[mfi]), - AMREX_REAL_ANYD(dx), - AMREX_MFITER_REDUCE_MIN(&dt)); - } - estdt_hydro = std::min(estdt_hydro, dt); - } + estdt_hydro = estdt_cfl(time); #ifdef RADIATION } @@ -1384,23 +1367,7 @@ Castro::estTimeStep (Real dt_old) if (diffuse_temp) { -#ifdef _OPENMP -#pragma omp parallel reduction(min:estdt_diffusion) -#endif - { - Real dt = max_dt / cfl; - - for (MFIter mfi(stateMF, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& box = mfi.tilebox(); - -#pragma gpu box(box) - ca_estdt_temp_diffusion(AMREX_INT_ANYD(box.loVect()), AMREX_INT_ANYD(box.hiVect()), - BL_TO_FORTRAN_ANYD(stateMF[mfi]), - AMREX_REAL_ANYD(dx), AMREX_MFITER_REDUCE_MIN(&dt)); - } - estdt_diffusion = std::min(estdt_diffusion, dt); - } + estdt_diffusion = estdt_temp_diffusion(); } ParallelDescriptor::ReduceRealMin(estdt_diffusion); diff --git a/Source/driver/Castro_F.H b/Source/driver/Castro_F.H index f15deb2cd8..51e36613cf 100644 --- a/Source/driver/Castro_F.H +++ b/Source/driver/Castro_F.H @@ -160,17 +160,6 @@ extern "C" (const int* lo, const int* hi, BL_FORT_FAB_ARG_3D(S_new), const int verbose); - void ca_estdt - (const int* lo, const int* hi, - const BL_FORT_FAB_ARG_3D(state), - const amrex::Real* dx, amrex::Real* dt); - -#ifdef DIFFUSION - void ca_estdt_temp_diffusion - (const int* lo, const int* hi, - const BL_FORT_FAB_ARG_3D(state), - const amrex::Real* dx, amrex::Real* dt); -#endif #ifdef RADIATION void ca_estdt_rad diff --git a/Source/driver/Make.package b/Source/driver/Make.package index 060b8b76a4..aad0e952a3 100644 --- a/Source/driver/Make.package +++ b/Source/driver/Make.package @@ -62,3 +62,4 @@ endif ca_F90EXE_sources += prob_params_nd.F90 ca_F90EXE_sources += Tagging_nd.F90 ca_F90EXE_sources += timestep_nd.F90 +CEXE_sources += timestep.cpp diff --git a/Source/driver/timestep.cpp b/Source/driver/timestep.cpp new file mode 100644 index 0000000000..6bc8eff3c6 --- /dev/null +++ b/Source/driver/timestep.cpp @@ -0,0 +1,225 @@ +#include "Castro.H" +#include "Castro_F.H" + +#ifdef DIFFUSION +#include "conductivity.H" +#endif + +using namespace amrex; + +Real +Castro::estdt_cfl(const Real time) +{ + + // Courant-condition limited timestep + + GpuArray center; + ca_get_center(center.begin()); + +#ifdef ROTATION + GpuArray omega; + get_omega(time, omega.begin()); +#endif + + const auto dx = geom.CellSizeArray(); + + ReduceOps reduce_op; + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + + const MultiFab& stateMF = get_new_data(State_Type); + + const int ltime_integration_method = time_integration_method; +#ifdef ROTATION + const int ldo_rotation = do_rotation; + const int lstate_in_rotating_frame = state_in_rotating_frame; +#endif + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(stateMF, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box& box = mfi.tilebox(); + + auto u = stateMF.array(mfi); + + reduce_op.eval(box, reduce_data, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept -> ReduceTuple + { + + Real rhoInv = 1.0_rt / u(i,j,k,URHO); + + eos_t eos_state; + eos_state.rho = u(i,j,k,URHO); + eos_state.T = u(i,j,k,UTEMP); + eos_state.e = u(i,j,k,UEINT) * rhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = u(i,j,k,UFS+n) * rhoInv; + } + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = u(i,j,k,UFX+n) * rhoInv; + } + + eos(eos_input_re, eos_state); + + // Compute velocity and then calculate CFL timestep. + + Real ux = u(i,j,k,UMX) * rhoInv; + Real uy = u(i,j,k,UMY) * rhoInv; + Real uz = u(i,j,k,UMZ) * rhoInv; + +#ifdef ROTATION + if (ldo_rotation == 1 && lstate_in_rotating_frame != 1) { + Real vel[3]; + vel[0] = ux; + vel[1] = uy; + vel[2] = uz; + + GeometryData geomdata = geom.data(); + + inertial_to_rotational_velocity_c(i, j, k, geomdata, + center.begin(), omega.begin(), time, vel); + + ux = vel[0]; + uy = vel[1]; + uz = vel[2]; + } +#endif + + Real c = eos_state.cs; + + Real dt1 = dx[0]/(c + std::abs(ux)); + + Real dt2; +#if AMREX_SPACEDIM >= 2 + dt2 = dx[1]/(c + std::abs(uy)); +#else + dt2 = dt1; +#endif + + Real dt3; +#if AMREX_SPACEDIM == 3 + dt3 = dx[2]/(c + std::abs(uz)); +#else + dt3 = dt1; +#endif + + // The CTU method has a less restrictive timestep than MOL-based + // schemes (including the true SDC). Since the simplified SDC + // solver is based on CTU, we can use its timestep. + if (ltime_integration_method == 0 || ltime_integration_method == 3) { + return {amrex::min(dt1, dt2, dt3)}; + + } else { + // method of lines-style constraint is tougher + Real dt_tmp = 1.0_rt/dt1; +#if AMREX_SPACEIM >= 2 + dt_tmp += 1.0_rt/dt2; +#endif +#if AMREX_SPACEDIM == 3 + dt_tmp += 1.0_rt/dt3; +#endif + + return 1.0_rt/dt_tmp; + } + + }); + + } + + ReduceTuple hv = reduce_data.value(); + Real estdt_hydro = amrex::get<0>(hv); + + return estdt_hydro; + +} + +#ifdef DIFFUSION + +Real +Castro::estdt_temp_diffusion(void) +{ + + // Diffusion-limited timestep + // + // dt < 0.5 dx**2 / D + // where D = k/(rho c_v), and k is the conductivity + + const auto dx = geom.CellSizeArray(); + + ReduceOps reduce_op; + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + + const MultiFab& stateMF = get_new_data(State_Type); + + const Real ldiffuse_cutoff_density = diffuse_cutoff_density; + const Real lmax_dt = max_dt; + const Real lcfl = cfl; + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(stateMF, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box& box = mfi.tilebox(); + + auto ustate = stateMF.array(mfi); + + reduce_op.eval(box, reduce_data, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept -> ReduceTuple + { + + if (ustate(i,j,k,URHO) > ldiffuse_cutoff_density) { + + Real rho_inv = 1.0_rt/ustate(i,j,k,URHO); + + // we need cv + eos_t eos_state; + eos_state.rho = ustate(i,j,k,URHO); + eos_state.T = ustate(i,j,k,UTEMP); + eos_state.e = ustate(i,j,k,UEINT) * rho_inv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = ustate(i,j,k,UFS+n) * rho_inv; + } + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = ustate(i,j,k,UFX+n) * rho_inv; + } + + eos(eos_input_re, eos_state); + + // we also need the conductivity + conductivity(eos_state); + + // maybe we should check (and take action) on negative cv here? + Real D = eos_state.conductivity * rho_inv / eos_state.cv; + + Real dt1 = 0.5_rt * dx[0]*dx[0] / D; + + Real dt2; +#if AMREX_SPACEDIM >= 2 + dt2 = 0.5_rt * dx[1]*dx[1] / D; +#else + dt2 = dt1; +#endif + + Real dt3; +#if AMREX_SPACEDIM >= 3 + dt3 = 0.5_rt * dx[2]*dx[2] / D; +#else + dt3 = dt1; +#endif + + return {amrex::min(dt1, dt2, dt3)}; + + } else { + return lmax_dt/lcfl; + } + }); + } + + ReduceTuple hv = reduce_data.value(); + Real estdt_diff = amrex::get<0>(hv); + + return estdt_diff; +} +#endif diff --git a/Source/driver/timestep_nd.F90 b/Source/driver/timestep_nd.F90 index 0e87efecc6..2a9f6ab186 100644 --- a/Source/driver/timestep_nd.F90 +++ b/Source/driver/timestep_nd.F90 @@ -9,115 +9,6 @@ module timestep_module contains - subroutine ca_estdt(lo,hi,u,u_lo,u_hi,dx,dt) bind(C, name="ca_estdt") - ! Courant-condition limited timestep - ! - ! .. note:: - ! Binds to C function ``ca_estdt`` - - use network, only: nspec, naux - use meth_params_module, only: NVAR, URHO, UMX, UMY, UMZ, UEINT, UTEMP, UFS, UFX, time_integration_method - use eos_module, only: eos - use eos_type_module, only: eos_t, eos_input_re - use prob_params_module, only: dim - use amrex_constants_module, ONLY : ONE -#ifdef ROTATION - use meth_params_module, only: do_rotation, state_in_rotating_frame - use rotation_module, only: inertial_to_rotational_velocity - use amrinfo_module, only: amr_time -#endif - use amrex_fort_module, only : rt => amrex_real - use reduction_module, only: reduce_min - - implicit none - - integer, intent(in) :: lo(3), hi(3) - integer, intent(in) :: u_lo(3), u_hi(3) - real(rt), intent(in) :: u(u_lo(1):u_hi(1),u_lo(2):u_hi(2),u_lo(3):u_hi(3),NVAR) - real(rt), intent(in) :: dx(3) - real(rt), intent(inout) :: dt - - real(rt) :: rhoInv, ux, uy, uz, c, dt1, dt2, dt3, dt_tmp - integer :: i, j, k - - type (eos_t) :: eos_state - -#ifdef ROTATION - real(rt) :: vel(3) -#endif - - !$gpu - - ! Call EOS for the purpose of computing sound speed - - do k = lo(3), hi(3) - do j = lo(2), hi(2) - do i = lo(1), hi(1) - rhoInv = ONE / u(i,j,k,URHO) - - eos_state % rho = u(i,j,k,URHO ) - eos_state % T = u(i,j,k,UTEMP) - eos_state % e = u(i,j,k,UEINT) * rhoInv - eos_state % xn = u(i,j,k,UFS:UFS+nspec-1) * rhoInv - eos_state % aux = u(i,j,k,UFX:UFX+naux-1) * rhoInv - - call eos(eos_input_re, eos_state) - - ! Compute velocity and then calculate CFL timestep. - - ux = u(i,j,k,UMX) * rhoInv - uy = u(i,j,k,UMY) * rhoInv - uz = u(i,j,k,UMZ) * rhoInv - -#ifdef ROTATION - if (do_rotation == 1 .and. state_in_rotating_frame /= 1) then - vel = [ux, uy, uz] - call inertial_to_rotational_velocity([i, j, k], amr_time, vel) - ux = vel(1) - uy = vel(2) - uz = vel(3) - endif -#endif - - c = eos_state % cs - - dt1 = dx(1)/(c + abs(ux)) - if (dim >= 2) then - dt2 = dx(2)/(c + abs(uy)) - else - dt2 = dt1 - endif - if (dim == 3) then - dt3 = dx(3)/(c + abs(uz)) - else - dt3 = dt1 - endif - - ! The CTU method has a less restrictive timestep than - ! MOL-based schemes (including the true SDC). Since the - ! simplified SDC solver is based on CTU, we can use its - ! timestep. - if (time_integration_method == 0 .or. time_integration_method == 3) then - call reduce_min(dt, min(dt1,dt2,dt3)) - else - ! method of lines-style constraint is tougher - dt_tmp = ONE/dt1 - if (dim >= 2) then - dt_tmp = dt_tmp + ONE/dt2 - endif - if (dim == 3) then - dt_tmp = dt_tmp + ONE/dt3 - endif - - call reduce_min(dt, ONE/dt_tmp) - endif - - enddo - enddo - enddo - - end subroutine ca_estdt - #ifdef REACTIONS subroutine ca_estdt_burning(lo, hi, & @@ -262,94 +153,6 @@ end subroutine ca_estdt_burning #endif -#ifdef DIFFUSION - - subroutine ca_estdt_temp_diffusion(lo, hi, & - state, s_lo, s_hi, & - dx, dt) bind(C, name="ca_estdt_temp_diffusion") - ! Diffusion-limited timestep - ! - ! .. note:: - ! Binds to C function ``ca_estdt_temp_diffusion`` - - use network, only: nspec, naux - use eos_module, only: eos - use eos_type_module, only: eos_input_re, eos_t - use meth_params_module, only: NVAR, URHO, UEINT, UTEMP, UFS, UFX, & - diffuse_cutoff_density - use prob_params_module, only: dim - use amrex_constants_module, only : ONE, HALF - use conductivity_module, only: conductivity - use amrex_fort_module, only: rt => amrex_real - use reduction_module, only: reduce_min - - implicit none - - integer, intent(in ) :: lo(3), hi(3) - integer, intent(in ) :: s_lo(3), s_hi(3) - real(rt), intent(in ) :: state(s_lo(1):s_hi(1),s_lo(2):s_hi(2),s_lo(3):s_hi(3),NVAR) - real(rt), intent(in ) :: dx(3) - real(rt), intent(inout) :: dt - - real(rt) :: dt1, dt2, dt3, rho_inv - integer :: i, j, k - real(rt) :: D - - type (eos_t) :: eos_state - - !$gpu - - ! dt < 0.5 dx**2 / D - ! where D = k/(rho c_v), and k is the conductivity - - do k = lo(3), hi(3) - do j = lo(2), hi(2) - do i = lo(1), hi(1) - - if (state(i,j,k,URHO) > diffuse_cutoff_density) then - - rho_inv = ONE/state(i,j,k,URHO) - - ! we need cv - eos_state % rho = state(i,j,k,URHO ) - eos_state % T = state(i,j,k,UTEMP) - eos_state % e = state(i,j,k,UEINT) * rho_inv - - eos_state % xn = state(i,j,k,UFS:UFS+nspec-1) * rho_inv - eos_state % aux = state(i,j,k,UFX:UFX+naux-1) * rho_inv - - call eos(eos_input_re, eos_state) - - ! we also need the conductivity - call conductivity(eos_state) - - ! maybe we should check (and take action) on negative cv here? - D = eos_state % conductivity*rho_inv/eos_state%cv - - dt1 = HALF*dx(1)**2/D - - if (dim >= 2) then - dt2 = HALF*dx(2)**2/D - else - dt2 = dt1 - endif - - if (dim == 3) then - dt3 = HALF*dx(3)**2/D - else - dt3 = dt1 - endif - - call reduce_min(dt, min(dt1,dt2,dt3)) - - endif - - enddo - enddo - enddo - - end subroutine ca_estdt_temp_diffusion -#endif subroutine ca_check_timestep(lo, hi, & diff --git a/Source/sources/Castro_sources_F.H b/Source/sources/Castro_sources_F.H index 73254542a7..e0c043b789 100644 --- a/Source/sources/Castro_sources_F.H +++ b/Source/sources/Castro_sources_F.H @@ -21,14 +21,6 @@ extern "C" void update_sponge_params(const amrex::Real* time); #endif - void ca_thermo_src - (const int* lo, const int* hi, - const BL_FORT_FAB_ARG_3D(state_old), - const BL_FORT_FAB_ARG_3D(state_new), - BL_FORT_FAB_ARG_3D(thermo_src), - const amrex::Real* prob_lo, const amrex::Real* dx, - const amrex::Real time, const amrex::Real dt); - void ca_ext_src (const int* lo, const int* hi, const BL_FORT_FAB_ARG_3D(old_state), diff --git a/Source/sources/Castro_thermo.cpp b/Source/sources/Castro_thermo.cpp index 842251ef65..48f1325711 100644 --- a/Source/sources/Castro_thermo.cpp +++ b/Source/sources/Castro_thermo.cpp @@ -55,22 +55,103 @@ Castro::fill_thermo_source (Real time, Real dt, MultiFab& state_old, MultiFab& state_new, MultiFab& thermo_src) { + + // Compute thermodynamic sources for the internal energy equation. + // At the moment, this is only the -p div{U} term in the internal + // energy equation, and only for method-of-lines integration, + // including the new SDC method (the `-` is because it is on the RHS + // of the equation) + const Real* dx = geom.CellSize(); const Real* prob_lo = geom.ProbLo(); + auto coord = geom.Coord(); + + #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(thermo_src, TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (MFIter mfi(thermo_src, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box& bx = mfi.tilebox(); + + Array4 const old_state = state_old.array(mfi); + Array4 const new_state = state_new.array(mfi); + Array4 const src = thermo_src.array(mfi); + + + AMREX_PARALLEL_FOR_3D(bx, i, j, k, { - const Box& bx = mfi.tilebox(); + // radius for non-Cartesian + Real rp = prob_lo[0] + (static_cast(i) + 1.5_rt)*dx[0]; + Real rm = prob_lo[0] + (static_cast(i) - 0.5_rt)*dx[0]; + Real r = 0.5_rt*(rm + rp); + + // compute -div{U} + if (coord == 0) { + src(i,j,k,UEINT) = -0.25_rt*((old_state(i+1,j,k,UMX)/old_state(i+1,j,k,URHO) + + new_state(i+1,j,k,UMX)/new_state(i+1,j,k,URHO)) - + (old_state(i-1,j,k,UMX)/old_state(i-1,j,k,URHO) + + new_state(i-1,j,k,UMX)/new_state(i-1,j,k,URHO)))/dx[0]; + + } else if (coord == 1) { + // axisymmetric + src(i,j,k,UEINT) = -0.25_rt*(rp*(old_state(i+1,j,k,UMX)/old_state(i+1,j,k,URHO) + + new_state(i+1,j,k,UMX)/new_state(i+1,j,k,URHO)) - + rm*(old_state(i-1,j,k,UMX)/old_state(i-1,j,k,URHO) + + new_state(i-1,j,k,UMX)/new_state(i-1,j,k,URHO)))/(r*dx[0]); + + } else if (coord == 2) { + // spherical + src(i,j,k,UEINT) = -0.25_rt*(rp*rp*(old_state(i+1,j,k,UMX)/old_state(i+1,j,k,URHO) + + new_state(i+1,j,k,UMX)/new_state(i+1,j,k,URHO)) - + rm*rm*(old_state(i-1,j,k,UMX)/old_state(i-1,j,k,URHO) + + new_state(i-1,j,k,UMX)/new_state(i-1,j,k,URHO)))/(r*r*dx[0]); + } + +#if BL_SPACEDIM >= 2 + src(i,j,k,UEINT) += -0.25_rt*((old_state(i,j+1,k,UMY)/old_state(i,j+1,k,URHO) + + new_state(i,j+1,k,UMY)/new_state(i,j+1,k,URHO)) - + (old_state(i,j-1,k,UMY)/old_state(i,j-1,k,URHO) + + new_state(i,j-1,k,UMY)/new_state(i,j-1,k,URHO)))/dx[1]; +#endif +#if BL_SPACEDIM == 3 + src(i,j,k,UEINT) += -0.25_rt*((old_state(i,j,k+1,UMZ)/old_state(i,j,k+1,URHO) + + new_state(i,j,k+1,UMZ)/new_state(i,j,k+1,URHO)) - + (old_state(i,j,k-1,UMZ)/old_state(i,j,k-1,URHO) + + new_state(i,j,k-1,UMZ)/new_state(i,j,k-1,URHO)))/dx[2]; +#endif -#pragma gpu box(bx) - ca_thermo_src(AMREX_INT_ANYD(bx.loVect()), AMREX_INT_ANYD(bx.hiVect()), - BL_TO_FORTRAN_ANYD(state_old[mfi]), - BL_TO_FORTRAN_ANYD(state_new[mfi]), - BL_TO_FORTRAN_ANYD(thermo_src[mfi]), - AMREX_REAL_ANYD(prob_lo),AMREX_REAL_ANYD(dx),time,dt); - } + // we now need the pressure -- we will assume that the + // temperature is consistent with the input state + eos_t eos_state_old; + eos_state_old.rho = old_state(i,j,k,URHO); + eos_state_old.T = old_state(i,j,k,UTEMP); + for (int n = 0; n < NumSpec; n++) { + eos_state_old.xn[n] = old_state(i,j,k,UFS+n)/old_state(i,j,k,URHO); + } + for (int n = 0; n < NumAux; n++) { + eos_state_old.aux[n] = old_state(i,j,k,UFX+n)/old_state(i,j,k,URHO); + } + + eos(eos_input_rt, eos_state_old); + + eos_t eos_state_new; + eos_state_new.rho = new_state(i,j,k,URHO); + eos_state_new.T = new_state(i,j,k,UTEMP); + for (int n = 0; n < NumSpec; n++) { + eos_state_new.xn[n] = new_state(i,j,k,UFS+n)/new_state(i,j,k,URHO); + } + for (int n = 0; n < NumAux; n++) { + eos_state_new.aux[n] = new_state(i,j,k,UFX+n)/new_state(i,j,k,URHO); + } + + eos(eos_input_rt, eos_state_new); + + // final source term, -p div{U} + src(i,j,k,UEINT) = 0.5_rt*(eos_state_old.p + eos_state_new.p)*src(i,j,k,UEINT); + + }); + } } diff --git a/Source/sources/Make.package b/Source/sources/Make.package index 570d8f17b0..895373f21a 100644 --- a/Source/sources/Make.package +++ b/Source/sources/Make.package @@ -10,7 +10,6 @@ CEXE_sources += Castro_thermo.cpp ca_f90EXE_sources += update_sponge_params.f90 ca_F90EXE_sources += sponge_nd.F90 -ca_F90EXE_sources += thermo_nd.F90 ca_F90EXE_sources += ext_src_nd.F90 diff --git a/Source/sources/thermo_nd.F90 b/Source/sources/thermo_nd.F90 deleted file mode 100644 index 5fdc869771..0000000000 --- a/Source/sources/thermo_nd.F90 +++ /dev/null @@ -1,120 +0,0 @@ -module thermo_sources - ! Functions for implementing the source terms to the internal energy - ! equation, for those integration methods where we treat it as a - ! source and do not explicitly discretize it in the conservative - ! update. - - implicit none - -contains - - subroutine ca_thermo_src(lo, hi, & - old_state, os_lo, os_hi, & - new_state, ns_lo, ns_hi, & - src, src_lo, src_hi, problo, dx, time, dt) bind(C, name="ca_thermo_src") - ! Compute thermodynamic sources for the internal energy equation. - ! At the moment, this is only the -p div{U} term in the internal energy - ! equation, and only for method-of-lines integration, including the new - ! SDC method (the `-` is because it is on the RHS of the equation) - - use amrex_constants_module, only: ZERO, HALF, FOURTH - use meth_params_module, only : NVAR, URHO, UMX, UMY, UMZ, UTEMP, UFS, UFX, UEINT, NSRC - use amrex_fort_module, only : rt => amrex_real - use network, only : nspec, naux - use eos_type_module, only : eos_t, eos_input_rt - use eos_module, only : eos - use prob_params_module, only : coord_type - - implicit none - - integer, intent(in) :: lo(3),hi(3) ! bounds of the box to operate on - integer, intent(in) :: os_lo(3),os_hi(3) ! bounds of the old_state array - integer, intent(in) :: ns_lo(3),ns_hi(3) ! bounds of the new_state array - integer, intent(in) :: src_lo(3),src_hi(3) ! bounds of the src array - real(rt), intent(in) :: old_state(os_lo(1):os_hi(1),os_lo(2):os_hi(2),os_lo(3):os_hi(3),NVAR) ! the old-time hydrodynamic conserved state - real(rt), intent(in) :: new_state(ns_lo(1):ns_hi(1),ns_lo(2):ns_hi(2),ns_lo(3):ns_hi(3),NVAR) ! the new-time hydrodynamic conserved state - real(rt), intent(inout) :: src(src_lo(1):src_hi(1),src_lo(2):src_hi(2),src_lo(3):src_hi(3),NSRC) ! the source terms for the conserved variables - real(rt), intent(in) :: problo(3) ! physical coordinates of the lower left corner of the domain - real(rt), intent(in) :: dx(3) ! grid spacing - real(rt), intent(in), value :: time ! current simulation time - real(rt), intent(in), value :: dt ! current timestep - - integer :: i, j, k - - type(eos_t) :: eos_state_old, eos_state_new - real(rt) :: r, rp, rm - - !$gpu - - do k = lo(3), hi(3) - do j = lo(2), hi(2) - do i = lo(1), hi(1) - - ! radius for non-Cartesian - rp = problo(1) + (dble(i) + 1.5_rt)*dx(1) - rm = problo(1) + (dble(i) - HALF)*dx(1) - r = HALF*(rm + rp) - - ! compute -div{U} - if (coord_type == 0) then - src(i,j,k,UEINT) = -FOURTH*((old_state(i+1,j,k,UMX)/old_state(i+1,j,k,URHO) + & - new_state(i+1,j,k,UMX)/new_state(i+1,j,k,URHO)) - & - (old_state(i-1,j,k,UMX)/old_state(i-1,j,k,URHO) + & - new_state(i-1,j,k,UMX)/new_state(i-1,j,k,URHO)))/dx(1) - - else if (coord_type == 1) then - ! axisymmetric - src(i,j,k,UEINT) = -FOURTH*(rp*(old_state(i+1,j,k,UMX)/old_state(i+1,j,k,URHO) + & - new_state(i+1,j,k,UMX)/new_state(i+1,j,k,URHO)) - & - rm*(old_state(i-1,j,k,UMX)/old_state(i-1,j,k,URHO) + & - new_state(i-1,j,k,UMX)/new_state(i-1,j,k,URHO)))/(r*dx(1)) - - else if (coord_type == 2) then - ! spherical - src(i,j,k,UEINT) = -FOURTH*(rp**2*(old_state(i+1,j,k,UMX)/old_state(i+1,j,k,URHO) + & - new_state(i+1,j,k,UMX)/new_state(i+1,j,k,URHO)) - & - rm**2*(old_state(i-1,j,k,UMX)/old_state(i-1,j,k,URHO) + & - new_state(i-1,j,k,UMX)/new_state(i-1,j,k,URHO)))/(r**2*dx(1)) - endif - -#if BL_SPACEDIM >= 2 - src(i,j,k,UEINT) = src(i,j,k,UEINT) - & - FOURTH*((old_state(i,j+1,k,UMY)/old_state(i,j+1,k,URHO) + & - new_state(i,j+1,k,UMY)/new_state(i,j+1,k,URHO)) - & - (old_state(i,j-1,k,UMY)/old_state(i,j-1,k,URHO) + & - new_state(i,j-1,k,UMY)/new_state(i,j-1,k,URHO)))/dx(2) -#endif -#if BL_SPACEDIM == 3 - src(i,j,k,UEINT) = src(i,j,k,UEINT) - & - FOURTH*((old_state(i,j,k+1,UMZ)/old_state(i,j,k+1,URHO) + & - new_state(i,j,k+1,UMZ)/new_state(i,j,k+1,URHO)) - & - (old_state(i,j,k-1,UMZ)/old_state(i,j,k-1,URHO) + & - new_state(i,j,k-1,UMZ)/new_state(i,j,k-1,URHO)))/dx(3) -#endif - - ! we now need the pressure -- we will assume that the - ! temperature is consistent with the input state - eos_state_old % rho = old_state(i,j,k,URHO) - eos_state_old % T = old_state(i,j,k,UTEMP) - eos_state_old % xn(:) = old_state(i,j,k,UFS:UFS-1+nspec)/old_state(i,j,k,URHO) - eos_state_old % aux(:) = old_state(i,j,k,UFX:UFX+naux-1)/old_state(i,j,k,URHO) - - call eos(eos_input_rt, eos_state_old) - - eos_state_new % rho = new_state(i,j,k,URHO) - eos_state_new % T = new_state(i,j,k,UTEMP) - eos_state_new % xn(:) = new_state(i,j,k,UFS:UFS-1+nspec)/new_state(i,j,k,URHO) - eos_state_new % aux(:) = new_state(i,j,k,UFX:UFX+naux-1)/new_state(i,j,k,URHO) - - call eos(eos_input_rt, eos_state_new) - - ! final source term, -p div{U} - src(i,j,k,UEINT) = HALF*(eos_state_old % p + eos_state_new % p)*src(i,j,k,UEINT) - - enddo - enddo - enddo - - end subroutine ca_thermo_src - -end module thermo_sources