diff --git a/docs/Models/linear-euler-3d-model.md b/docs/Models/linear-euler-3d-model.md index 00b19d2cc..c7307bc97 100644 --- a/docs/Models/linear-euler-3d-model.md +++ b/docs/Models/linear-euler-3d-model.md @@ -20,11 +20,12 @@ $$ u \\ v \\ w \\ - p + p \\ + c \end{pmatrix} $$ -where $\rho$ is a density anomaly, referenced to the density $\rho_0$, $u$, $v$, and $w$ are the $x$, $y$, and $z$ components of the fluid velocity (respectively), and $p$ is the pressure. When we assume an ideal gas, and a motionless background state, the conservative fluxes are +where $\rho$ is a density anomaly, referenced to the density $\rho_0$, $u$, $v$, and $w$ are the $x$, $y$, and $z$ components of the fluid velocity (respectively), and $p$ is the pressure. The sound speed $c$ is carried as a solution variable so that heterogeneous (spatially varying) media are supported; it has zero flux and zero source, so it is static in time and is set by the initial condition (mirroring the 2-D model). When we assume an ideal gas, and a motionless background state, the conservative fluxes are $$ \overleftrightarrow{f} = @@ -33,11 +34,12 @@ $$ p \hat{x} \\ p \hat{y} \\ p \hat{z} \\ - \rho_0c^2(u \hat{x} + v \hat{y} + w \hat{z}) + \rho_0c^2(u \hat{x} + v \hat{y} + w \hat{z}) \\ + 0 \end{pmatrix} $$ -where $c$ is the (constant) speed of sound. The source term is set to zero. +The source term is set to zero. $$ \vec{q} = \vec{0} @@ -46,44 +48,35 @@ $$ To track stability of the Euler equation, the total entropy function is $$ - e = \frac{1}{2} \int_V u^2 + v^2 + \frac{p}{\rho_0 c^2} \hspace{1mm} dV + e = \frac{1}{2} \int_V \rho_0 (u^2 + v^2 + w^2) + \frac{p^2}{\rho_0 c^2} \hspace{1mm} dV $$ ## Implementation -The Linear Euler 3D model is implemented as a type extension of the [`DGModel3D` class](../ford/type/dgmodel3d_t.html). The [`LinearEuler3D_t` class](../ford/type/lineareuler3d_t.html) adds parameters for the reference density and the speed speed of sound and overrides the `SetMetadata`, `entropy_func`, `flux3d`, and `riemannflux3d` type-bound procedures. +The Linear Euler 3D model is implemented as a type extension of the [`DGModel3D` class](../ford/type/dgmodel3d_t.html). The [`LinearEuler3D_t` class](../ford/type/lineareuler3d_t.html) adds parameters for the reference density and the reference speed of sound and overrides the `SetMetadata`, `entropy_func`, `flux3d`, and `riemannflux3d` type-bound procedures. ### Riemann Solver -The `LinearEuler3D` class is defined using the conservative form of the conservation law. The Riemman solver for the hyperbolic part of Euler equation is the local Lax Friedrichs upwind riemann solver +The `LinearEuler3D` class is defined using the conservative form of the conservation law. The Riemann solver for the hyperbolic part of the Euler equation is the impedance-matched (characteristic/Godunov) upwind flux, identical in form to the 2-D model's. With the acoustic impedances $Z_L = \rho_0 c_L$ and $Z_R = \rho_0 c_R$ evaluated from the per-node sound speed on either side of the face, the interface normal velocity and pressure are $$ - \overleftrightarrow{f}_h^* \cdot \hat{n} = \frac{1}{2}( \overleftrightarrow{f}_L \cdot \hat{n} + \overleftrightarrow{f}_R \cdot \hat{n} + c (\vec{s}_L - \vec{s}_R)) + u_n^* = \frac{Z_L u_{n,L} + Z_R u_{n,R} + (p_L - p_R)}{Z_L + Z_R}, \qquad + p^* = \frac{Z_R p_L + Z_L p_R + Z_L Z_R (u_{n,L} - u_{n,R})}{Z_L + Z_R} $$ -where +and the normal flux is $$ - \overleftrightarrow{f}_L \cdot \hat{n} = + \overleftrightarrow{f}_h^* \cdot \hat{n} = \begin{pmatrix} - \rho_0(u_L n_x + v_L n_y + w_L n_z) \\ - p_L n_x \\ - p_L n_y \\ - p_L n_z \\ - \rho_0c^2(u_L n_x + v_L n_y + w_L n_z) - \end{pmatrix} -$$ - -$$ - \overleftrightarrow{f}_R \cdot \hat{n} = - \begin{pmatrix} - \rho_0(u_R n_x + v_R n_y + w_R n_z) \\ - p_R n_x \\ - p_R n_y \\ - p_R n_z \\ - \rho_0c^2(u_R n_x + v_R n_y + w_R n_z) - \end{pmatrix} + \rho_0 u_n^* \\ + p^* n_x / \rho_0 \\ + p^* n_y / \rho_0 \\ + p^* n_z / \rho_0 \\ + \rho_0 \overline{c^2} u_n^* \\ + 0 + \end{pmatrix}, \qquad \overline{c^2} = \frac{1}{2}(c_L^2 + c_R^2) $$ -The details for this implementation can be found in [`self_lineareuler3d_t.f90`](../ford/sourcefile/self_lineareuler3d_t.f90.html) +Because the interface states are resolved with the one-sided impedances, material interfaces in a heterogeneous sound-speed field produce the physically correct transmission and reflection. The details for this implementation can be found in [`self_lineareuler3d_t.f90`](../ford/sourcefile/self_lineareuler3d_t.f90.html) ### Boundary conditions When initializing the mesh for your Euler 3D equation solver, you can change the boundary conditions to diff --git a/src/SELF_LinearEuler3D_t.f90 b/src/SELF_LinearEuler3D_t.f90 index e51ca033d..473875637 100644 --- a/src/SELF_LinearEuler3D_t.f90 +++ b/src/SELF_LinearEuler3D_t.f90 @@ -27,17 +27,20 @@ module self_LinearEuler3D_t !! This module defines a class that can be used to solve the Linear Euler !! equations in 3-D. The Linear Euler Equations, here, are the Euler equations -!! linearized about a motionless background state. +!! linearized about a motionless background state. The sound speed is carried +!! as a solution variable (static in time, possibly spatially varying) so that +!! heterogeneous media are supported, mirroring the 2-D model. +!! +!! The solution variables are !! -!! The conserved variables are - !! \begin{equation} !! \vec{s} = \begin{pmatrix} !! \rho \\ !! u \\ !! v \\ !! w \\ -!! p +!! p \\ +!! c !! \end{pmatrix} !! \end{equation} !! @@ -49,31 +52,36 @@ module self_LinearEuler3D_t !! \frac{p}{\rho_0} \hat{x} \\ !! \frac{p}{\rho_0} \hat{y} \\ !! \frac{p}{\rho_0} \hat{z} \\ -!! c^2 \rho_0 ( u \hat{x} + v \hat{y} + w \hat{z}) +!! c^2 \rho_0 ( u \hat{x} + v \hat{y} + w \hat{z}) \\ +!! 0 !! \end{pmatrix} !! \end{equation} !! -!! and the source terms are null. +!! and the source terms are null. The sound speed variable has zero flux and +!! zero source; it is set by the initial condition and preserved in time. +!! The Riemann flux is the impedance-matched (characteristic) flux, identical +!! in form to the 2-D model's; see riemannflux3d_LinearEuler3D_t. !! use self_model use self_dgmodel3D use self_mesh + use SELF_BoundaryConditions implicit none type,extends(dgmodel3D) :: LinearEuler3D_t ! Add any additional attributes here that are specific to your model real(prec) :: rho0 = 1.0_prec ! Reference density - real(prec) :: c = 1.0_prec ! Sound speed + real(prec) :: c = 1.0_prec ! Reference sound speed (used to fill variable 6 in initial conditions) real(prec) :: g = 0.0_prec ! gravitational acceleration (y-direction only) contains procedure :: SourceMethod => sourcemethod_LinearEuler3D_t procedure :: SetNumberOfVariables => SetNumberOfVariables_LinearEuler3D_t procedure :: SetMetadata => SetMetadata_LinearEuler3D_t + procedure :: AdditionalInit => AdditionalInit_LinearEuler3D_t procedure :: entropy_func => entropy_func_LinearEuler3D_t - !procedure :: hbc3D_NoNormalFlow => hbc3D_NoNormalFlow_LinearEuler3D_t procedure :: flux3D => flux3D_LinearEuler3D_t procedure :: riemannflux3D => riemannflux3D_LinearEuler3D_t procedure :: SphericalSoundWave => SphericalSoundWave_LinearEuler3D_t @@ -86,7 +94,7 @@ subroutine SetNumberOfVariables_LinearEuler3D_t(this) implicit none class(LinearEuler3D_t),intent(inout) :: this - this%nvar = 5 + this%nvar = 6 endsubroutine SetNumberOfVariables_LinearEuler3D_t @@ -109,38 +117,65 @@ subroutine SetMetadata_LinearEuler3D_t(this) call this%solution%SetName(5,"P") ! Pressure call this%solution%SetUnits(5,"kg⋅m⁻¹⋅s⁻²") + call this%solution%SetName(6,"c") ! Sound speed (static; possibly heterogeneous) + call this%solution%SetUnits(6,"m⋅s⁻¹") + endsubroutine SetMetadata_LinearEuler3D_t + subroutine AdditionalInit_LinearEuler3D_t(this) + !! Register the (CPU) radiation boundary condition. GPU builds overwrite + !! this registration with the device kernel in AdditionalInit_LinearEuler3D. + implicit none + class(LinearEuler3D_t),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + bcfunc => hbc3d_Radiation_LinearEuler3D + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_RADIATION,"radiation",bcfunc) + + endsubroutine AdditionalInit_LinearEuler3D_t + + subroutine hbc3d_Radiation_LinearEuler3D(bc,mymodel) + !! Radiation BC: zero acoustic perturbation in the exterior state; the + !! sound speed (variable 6) is copied from the interior side so the + !! Riemann solver sees a consistent c. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,j,iEl,s + + select type(m => mymodel) + class is(LinearEuler3D_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + do j = 1,m%solution%interp%N+1 + do i = 1,m%solution%interp%N+1 + m%solution%extBoundary(i,j,s,iEl,1:5) = 0.0_prec + m%solution%extBoundary(i,j,s,iEl,6) = m%solution%boundary(i,j,s,iEl,6) ! c preserved + enddo + enddo + enddo + endselect + + endsubroutine hbc3d_Radiation_LinearEuler3D + pure function entropy_func_LinearEuler3D_t(this,s) result(e) !! The entropy function is the sum of kinetic and internal energy !! For the linear model, this is !! !! \begin{equation} - !! e = \frac{1}{2} \left( \rho_0*( u^2 + v^2 ) + \frac{P^2}{\rho_0 c^2} \right) + !! e = \frac{1}{2} \left( \rho_0*( u^2 + v^2 + w^2 ) + \frac{P^2}{\rho_0 c^2} \right) class(LinearEuler3D_t),intent(in) :: this real(prec),intent(in) :: s(1:this%nvar) real(prec) :: e - e = 0.5_prec*this%rho0*(s(2)*s(2)+s(3)*(3)+s(4)*s(4))+ & - 0.5_prec*(s(5)*s(5)/(this%rho0*this%c*this%c)) + e = 0.5_prec*this%rho0*(s(2)*s(2)+s(3)*s(3)+s(4)*s(4))+ & + 0.5_prec*(s(5)*s(5)/(this%rho0*s(6)*s(6))) endfunction entropy_func_LinearEuler3D_t - ! pure function hbc3D_NoNormalFlow_LinearEuler3D_t(this,s,nhat) result(exts) - ! class(LinearEuler3D_t),intent(in) :: this - ! real(prec),intent(in) :: s(1:this%nvar) - ! real(prec),intent(in) :: nhat(1:2) - ! real(prec) :: exts(1:this%nvar) - ! ! Local - ! integer :: ivar - - ! exts(1) = s(1) ! density - ! exts(2) = (nhat(2)**2-nhat(1)**2)*s(2)-2.0_prec*nhat(1)*nhat(2)*s(3) ! u - ! exts(3) = (nhat(1)**2-nhat(2)**2)*s(3)-2.0_prec*nhat(1)*nhat(2)*s(2) ! v - ! exts(4) = (nhat(1)**2-nhat(2)**2)*s(3)-2.0_prec*nhat(1)*nhat(2)*s(2) ! w - ! exts(5) = s(4) ! p - - ! endfunction hbc3D_NoNormalFlow_LinearEuler3D_t subroutine sourcemethod_LinearEuler3D_t(this) implicit none class(LinearEuler3D_t),intent(inout) :: this @@ -171,16 +206,28 @@ pure function flux3D_LinearEuler3D_t(this,s,dsdx) result(flux) flux(4,2) = 0.0_prec ! z-velocity, y flux; 0 flux(4,3) = s(5)/this%rho0 ! z-velocity, z flux; p/rho0 - flux(5,1) = this%c*this%c*this%rho0*s(2) ! pressure, x flux : rho0*c^2*u - flux(5,2) = this%c*this%c*this%rho0*s(3) ! pressure, y flux : rho0*c^2*v - flux(5,3) = this%c*this%c*this%rho0*s(4) ! pressure, y flux : rho0*c^2*w + flux(5,1) = s(6)*s(6)*this%rho0*s(2) ! pressure, x flux : rho0*c^2*u + flux(5,2) = s(6)*s(6)*this%rho0*s(3) ! pressure, y flux : rho0*c^2*v + flux(5,3) = s(6)*s(6)*this%rho0*s(4) ! pressure, z flux : rho0*c^2*w + + flux(6,1) = 0.0_prec ! sound speed, x flux; 0 (c held fixed in time) + flux(6,2) = 0.0_prec ! sound speed, y flux; 0 + flux(6,3) = 0.0_prec ! sound speed, z flux; 0 if(.false.) flux(1,1) = flux(1,1)+dsdx(1,1) ! suppress unused-dummy-argument warning endfunction flux3D_LinearEuler3D_t pure function riemannflux3D_LinearEuler3D_t(this,sL,sR,dsdx,nhat) result(flux) - !! Uses a local lax-friedrich's upwind flux - !! The max eigenvalue is taken as the sound speed + !! Impedance-matched (characteristic/Godunov) Riemann flux, identical in + !! form to the 2-D model's. The interface states are resolved with the + !! acoustic impedances Z = rho0*c on either side, so material interfaces + !! in a heterogeneous sound-speed field are handled with the physically + !! correct transmission/reflection: + !! + !! un* = (ZL*unL + ZR*unR + (pL - pR)) / (ZL + ZR) + !! p* = (ZR*pL + ZL*pR + ZL*ZR*(unL - unR)) / (ZL + ZR) + !! + !! The sound speed variable carries zero flux. class(LinearEuler3D_t),intent(in) :: this real(prec),intent(in) :: sL(1:this%nvar) real(prec),intent(in) :: sR(1:this%nvar) @@ -188,33 +235,29 @@ pure function riemannflux3D_LinearEuler3D_t(this,sL,sR,dsdx,nhat) result(flux) real(prec),intent(in) :: nhat(1:3) real(prec) :: flux(1:this%nvar) ! Local - real(prec) :: fL(1:this%nvar) - real(prec) :: fR(1:this%nvar) - real(prec) :: u,v,w,p,c,rho0 - - u = sL(2) - v = sL(3) - w = sL(4) - p = sL(5) + real(prec) :: rho0,cL,cR,ZL,ZR,unL,unR,pL,pR,un_star,p_star,c2_avg + rho0 = this%rho0 - c = this%c - fL(1) = rho0*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! density - fL(2) = p*nhat(1)/rho0 ! u - fL(3) = p*nhat(2)/rho0 ! v - fL(4) = p*nhat(3)/rho0 ! w - fL(5) = rho0*c*c*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! pressure - - u = sR(2) - v = sR(3) - w = sR(4) - p = sR(5) - fR(1) = rho0*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! density - fR(2) = p*nhat(1)/rho0 ! u - fR(3) = p*nhat(2)/rho0 ! v' - fR(4) = p*nhat(3)/rho0 ! w - fR(5) = rho0*c*c*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! pressure - - flux(1:5) = 0.5_prec*(fL(1:5)+fR(1:5))+c*(sL(1:5)-sR(1:5)) + cL = sL(6) + cR = sR(6) + ZL = rho0*cL + ZR = rho0*cR + + unL = sL(2)*nhat(1)+sL(3)*nhat(2)+sL(4)*nhat(3) + unR = sR(2)*nhat(1)+sR(3)*nhat(2)+sR(4)*nhat(3) + pL = sL(5) + pR = sR(5) + + un_star = (ZL*unL+ZR*unR+(pL-pR))/(ZL+ZR) + p_star = (ZR*pL+ZL*pR+ZL*ZR*(unL-unR))/(ZL+ZR) + c2_avg = 0.5_prec*(cL*cL+cR*cR) + + flux(1) = rho0*un_star + flux(2) = p_star*nhat(1)/rho0 + flux(3) = p_star*nhat(2)/rho0 + flux(4) = p_star*nhat(3)/rho0 + flux(5) = rho0*c2_avg*un_star + flux(6) = 0.0_prec if(.false.) flux(1) = flux(1)+dsdx(1,1) ! suppress unused-dummy-argument warning endfunction riemannflux3D_LinearEuler3D_t @@ -260,6 +303,7 @@ subroutine SphericalSoundWave_LinearEuler3D_t(this,rhoprime,Lr,x0,y0,z0) this%solution%interior(i,j,k,iEl,3) = 0.0_prec this%solution%interior(i,j,k,iEl,4) = 0.0_prec this%solution%interior(i,j,k,iEl,5) = rho*this%c*this%c + this%solution%interior(i,j,k,iEl,6) = this%c ! uniform background sound speed enddo diff --git a/src/gpu/SELF_LinearEuler3D.cpp b/src/gpu/SELF_LinearEuler3D.cpp index 0d61687e0..36b51283d 100644 --- a/src/gpu/SELF_LinearEuler3D.cpp +++ b/src/gpu/SELF_LinearEuler3D.cpp @@ -29,46 +29,50 @@ #include "SELF_GPU_Macros.h" -__global__ void boundaryflux_LinearEuler3D_kernel(real *fb, real *extfb, real *nhat, real *nmag, real *flux, real rho0, real c, int ndof){ +// Impedance-matched (characteristic/Godunov) Riemann flux with a per-node +// sound speed carried as solution variable 6 (0-based index 5), mirroring the +// 2-D model. The interface states are resolved with the acoustic impedances +// Z = rho0*c on either side: +// un* = (ZL*unL + ZR*unR + (pL - pR)) / (ZL + ZR) +// p* = (ZR*pL + ZL*pR + ZL*ZR*(unL - unR)) / (ZL + ZR) +// Variable layout (0-based): 0=rho, 1=u, 2=v, 3=w, 4=p, 5=c. The sound-speed +// variable carries zero flux. +__global__ void boundaryflux_LinearEuler3D_kernel(real *fb, real *extfb, real *nhat, real *nmag, real *flux, real rho0, int ndof){ uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; - + if( idof < ndof ){ - real fl[5]; real nx = nhat[idof]; real ny = nhat[idof+ndof]; real nz = nhat[idof+2*ndof]; - real un = fb[idof + ndof]*nx + fb[idof + 2*ndof]*ny+fb[idof + 3*ndof]*nz; - real p = fb[idof + 4*ndof]; - - fl[0] = rho0*un; // density flux - fl[1] = p*nx/rho0; // x-momentum flux - fl[2] = p*ny/rho0; // y-momentum flux - fl[3] = p*nz/rho0; // z-momentum flux - fl[4] = rho0*c*c*un; // pressure flux - - real fr[5]; - un = extfb[idof + ndof]*nx + extfb[idof + 2*ndof]*ny+extfb[idof + 3*ndof]*nz; - p = extfb[idof + 4*ndof]; - - fr[0] = rho0*un; // density flux - fr[1] = p*nx/rho0; // x-momentum flux - fr[2] = p*ny/rho0; // y-momentum flux - fr[3] = p*nz/rho0; // y-momentum flux - fr[4] = rho0*c*c*un; // pressure flux - real nm = nmag[idof]; - flux[idof] = (0.5*(fl[0]+fr[0])+c*(fb[idof]-extfb[idof]))*nm; // density - flux[idof+ndof] = (0.5*(fl[1]+fr[1])+c*(fb[idof+ndof]-extfb[idof+ndof]))*nm; // u - flux[idof+2*ndof] = (0.5*(fl[2]+fr[2])+c*(fb[idof+2*ndof]-extfb[idof+2*ndof]))*nm; // v - flux[idof+3*ndof] = (0.5*(fl[3]+fr[3])+c*(fb[idof+3*ndof]-extfb[idof+3*ndof]))*nm; // w - flux[idof+4*ndof] = (0.5*(fl[4]+fr[4])+c*(fb[idof+4*ndof]-extfb[idof+4*ndof]))*nm; // p + + real cL = fb[idof + 5*ndof]; + real cR = extfb[idof + 5*ndof]; + real ZL = rho0*cL; + real ZR = rho0*cR; + + real unL = fb[idof + ndof]*nx + fb[idof + 2*ndof]*ny + fb[idof + 3*ndof]*nz; + real unR = extfb[idof + ndof]*nx + extfb[idof + 2*ndof]*ny + extfb[idof + 3*ndof]*nz; + real pL = fb[idof + 4*ndof]; + real pR = extfb[idof + 4*ndof]; + + real un_star = (ZL*unL + ZR*unR + (pL - pR))/(ZL + ZR); + real p_star = (ZR*pL + ZL*pR + ZL*ZR*(unL - unR))/(ZL + ZR); + real c2_avg = 0.5*(cL*cL + cR*cR); + + flux[idof] = (rho0*un_star)*nm; // density + flux[idof+ndof] = (p_star*nx/rho0)*nm; // u + flux[idof+2*ndof] = (p_star*ny/rho0)*nm; // v + flux[idof+3*ndof] = (p_star*nz/rho0)*nm; // w + flux[idof+4*ndof] = (rho0*c2_avg*un_star)*nm; // pressure + flux[idof+5*ndof] = 0.0; // sound speed (static) } } extern "C" { - void boundaryflux_LinearEuler3D_gpu(real *fb, real *extfb,real *nhat, real *nmag, real *flux, real rho0, real c, int N, int nel){ + void boundaryflux_LinearEuler3D_gpu(real *fb, real *extfb,real *nhat, real *nmag, real *flux, real rho0, int N, int nel, int nvar){ int threads_per_block = 256; uint32_t ndof = (N+1)*(N+1)*6*nel; int nblocks_x = ndof/threads_per_block +1; @@ -76,11 +80,11 @@ extern "C" dim3 nblocks(nblocks_x,1,1); dim3 nthreads(threads_per_block,1,1); - boundaryflux_LinearEuler3D_kernel<<>>(fb,extfb,nhat,nmag,flux,rho0,c,ndof); + boundaryflux_LinearEuler3D_kernel<<>>(fb,extfb,nhat,nmag,flux,rho0,ndof); } } - __global__ void fluxmethod_LinearEuler3D_gpukernel(real *solution, real *flux, real rho0, real c, int ndof, int nvar){ + __global__ void fluxmethod_LinearEuler3D_gpukernel(real *solution, real *flux, real rho0, int ndof, int nvar){ uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; if( idof < ndof ){ @@ -88,6 +92,7 @@ extern "C" real v = solution[idof + 2*ndof]; real w = solution[idof + 3*ndof]; real p = solution[idof + 4*ndof]; + real c = solution[idof + 5*ndof]; flux[idof + ndof*(0 + nvar*0)] = rho0*u; // density, x flux ; rho0*u flux[idof + ndof*(0 + nvar*1)] = rho0*v; // density, y flux ; rho0*v @@ -109,23 +114,29 @@ extern "C" flux[idof + ndof*(4 + nvar*1)] = c*c*rho0*v; // pressure, y flux : rho0*c^2*v flux[idof + ndof*(4 + nvar*2)] = c*c*rho0*w; // pressure, z flux : rho0*c^2*w + flux[idof + ndof*(5 + nvar*0)] = 0.0; // sound speed, x flux; 0 (c held fixed in time) + flux[idof + ndof*(5 + nvar*1)] = 0.0; // sound speed, y flux; 0 + flux[idof + ndof*(5 + nvar*2)] = 0.0; // sound speed, z flux; 0 + } } extern "C" { - void fluxmethod_LinearEuler3D_gpu(real *solution, real *flux, real rho0, real c, int N, int nel, int nvar){ + void fluxmethod_LinearEuler3D_gpu(real *solution, real *flux, real rho0, int N, int nel, int nvar){ int ndof = (N+1)*(N+1)*(N+1)*nel; int threads_per_block = 256; int nblocks_x = ndof/threads_per_block +1; - fluxmethod_LinearEuler3D_gpukernel<<>>(solution,flux,rho0,c,ndof,nvar); + fluxmethod_LinearEuler3D_gpukernel<<>>(solution,flux,rho0,ndof,nvar); } } // Radiation BC kernel for 3D Linear Euler -// Sets extBoundary = 0 on pre-filtered boundary faces +// Zeroes the acoustic perturbation (rho,u,v,w,p) in extBoundary on +// pre-filtered boundary faces; the sound speed (variable 5, 0-based) is +// copied from the interior side so face Riemann fluxes see a consistent c. __global__ void hbc3d_radiation_lineareuler3d_kernel( - real *extBoundary, + real *extBoundary, real *boundary, int *elements, int *sides, int nBoundaries, int N, int nel) { @@ -145,13 +156,14 @@ __global__ void hbc3d_radiation_lineareuler3d_kernel( extBoundary[SCB_3D_INDEX(i,j,s1,e1,2,N,nel)] = 0.0; extBoundary[SCB_3D_INDEX(i,j,s1,e1,3,N,nel)] = 0.0; extBoundary[SCB_3D_INDEX(i,j,s1,e1,4,N,nel)] = 0.0; + extBoundary[SCB_3D_INDEX(i,j,s1,e1,5,N,nel)] = boundary[SCB_3D_INDEX(i,j,s1,e1,5,N,nel)]; // c preserved } } extern "C" { void hbc3d_radiation_lineareuler3d_gpu( - real *extBoundary, + real *extBoundary, real *boundary, int *elements, int *sides, int nBoundaries, int N, int nel) { @@ -160,7 +172,7 @@ extern "C" int total_dofs = nBoundaries * dofs_per_face; int nblocks_x = total_dofs/threads_per_block + 1; hbc3d_radiation_lineareuler3d_kernel<<>>(extBoundary, + dim3(threads_per_block,1,1), 0, 0>>>(extBoundary, boundary, elements, sides, nBoundaries, N, nel); } } diff --git a/src/gpu/SELF_LinearEuler3D.f90 b/src/gpu/SELF_LinearEuler3D.f90 index 4da12fb41..aeec2a581 100644 --- a/src/gpu/SELF_LinearEuler3D.f90 +++ b/src/gpu/SELF_LinearEuler3D.f90 @@ -41,34 +41,34 @@ module self_LinearEuler3D endtype LinearEuler3D interface - subroutine hbc3d_radiation_lineareuler3d_gpu(extboundary, & + subroutine hbc3d_radiation_lineareuler3d_gpu(extboundary,boundary, & elements,sides,nBoundaries,N,nel) & bind(c,name="hbc3d_radiation_lineareuler3d_gpu") use iso_c_binding - type(c_ptr),value :: extboundary,elements,sides + type(c_ptr),value :: extboundary,boundary,elements,sides integer(c_int),value :: nBoundaries,N,nel endsubroutine hbc3d_radiation_lineareuler3d_gpu endinterface interface - subroutine fluxmethod_LinearEuler3D_gpu(solution,flux,rho0,c,N,nel,nvar) & + subroutine fluxmethod_LinearEuler3D_gpu(solution,flux,rho0,N,nel,nvar) & bind(c,name="fluxmethod_LinearEuler3D_gpu") use iso_c_binding use SELF_Constants type(c_ptr),value :: solution,flux - real(c_prec),value :: rho0,c + real(c_prec),value :: rho0 integer(c_int),value :: N,nel,nvar endsubroutine fluxmethod_LinearEuler3D_gpu endinterface interface - subroutine boundaryflux_LinearEuler3D_gpu(fb,fextb,nhat,nscale,flux,rho0,c,N,nel) & + subroutine boundaryflux_LinearEuler3D_gpu(fb,fextb,nhat,nscale,flux,rho0,N,nel,nvar) & bind(c,name="boundaryflux_LinearEuler3D_gpu") use iso_c_binding use SELF_Constants type(c_ptr),value :: fb,fextb,flux,nhat,nscale - real(c_prec),value :: rho0,c - integer(c_int),value :: N,nel + real(c_prec),value :: rho0 + integer(c_int),value :: N,nel,nvar endsubroutine boundaryflux_LinearEuler3D_gpu endinterface @@ -80,6 +80,10 @@ subroutine AdditionalInit_LinearEuler3D(this) ! Local procedure(SELF_bcMethod),pointer :: bcfunc + ! Register GPU-accelerated BC methods, overwriting the CPU versions + ! from the parent _t AdditionalInit + call AdditionalInit_LinearEuler3D_t(this) + bcfunc => hbc3d_Radiation_LinearEuler3D_GPU_wrapper call this%hyperbolicBCs%RegisterBoundaryCondition( & SELF_BC_RADIATION,"radiation",bcfunc) @@ -95,6 +99,7 @@ subroutine hbc3d_Radiation_LinearEuler3D_GPU_wrapper(bc,mymodel) if(bc%nBoundaries > 0) then call hbc3d_radiation_lineareuler3d_gpu( & m%solution%extBoundary_gpu, & + m%solution%boundary_gpu, & bc%elements_gpu,bc%sides_gpu, & bc%nBoundaries,m%solution%interp%N,m%solution%nElem) endif @@ -111,8 +116,8 @@ subroutine boundaryflux_LinearEuler3D(this) this%geometry%nhat%boundary_gpu, & this%geometry%nscale%boundary_gpu, & this%flux%boundarynormal_gpu, & - this%rho0,this%c,this%solution%interp%N, & - this%solution%nelem) + this%rho0,this%solution%interp%N, & + this%solution%nelem,this%solution%nvar) endsubroutine boundaryflux_LinearEuler3D @@ -122,7 +127,7 @@ subroutine fluxmethod_LinearEuler3D(this) call fluxmethod_LinearEuler3D_gpu(this%solution%interior_gpu, & this%flux%interior_gpu, & - this%rho0,this%c,this%solution%interp%N,this%solution%nelem, & + this%rho0,this%solution%interp%N,this%solution%nelem, & this%solution%nvar) endsubroutine fluxmethod_LinearEuler3D