From 9a47877b75a3b456e3f1ffd8402d52a741b097f8 Mon Sep 17 00:00:00 2001 From: Carlos Munoz Moncayo Date: Thu, 4 Jul 2024 12:23:29 +0300 Subject: [PATCH 1/7] Move source to examples directory --- examples/euler_gravity_3d/Makefile | 5 +- examples/euler_gravity_3d/__init__.py | 0 .../euler_gravity_3d/euler3d_mappedGrid.f90 | 679 ++++++++++++++++++ .../euler_gravity_3d/rising_hot_sphere.py | 12 +- .../rising_hot_sphere_spherical.py | 10 +- .../euler_gravity_3d/rpn3_euler_mapgrid.f90 | 312 ++++++++ .../euler_gravity_3d/rpt3_euler_mapgrid.f90 | 295 ++++++++ .../euler_gravity_3d/rptt3_euler_mapgrid.f90 | 306 ++++++++ 8 files changed, 1607 insertions(+), 12 deletions(-) create mode 100644 examples/euler_gravity_3d/__init__.py create mode 100644 examples/euler_gravity_3d/euler3d_mappedGrid.f90 create mode 100644 examples/euler_gravity_3d/rpn3_euler_mapgrid.f90 create mode 100644 examples/euler_gravity_3d/rpt3_euler_mapgrid.f90 create mode 100644 examples/euler_gravity_3d/rptt3_euler_mapgrid.f90 diff --git a/examples/euler_gravity_3d/Makefile b/examples/euler_gravity_3d/Makefile index 182410491..1ca2c7a93 100644 --- a/examples/euler_gravity_3d/Makefile +++ b/examples/euler_gravity_3d/Makefile @@ -1,15 +1,14 @@ F2PY = f2py -RIEMANN_SOURCE = ../../../riemann/src/ all: $(MAKE) euler_3d_gmap.so $(MAKE) mappedGrid.so -euler_3d_gmap.so: $(RIEMANN_SOURCE)/rpn3_euler_mapgrid.f90 $(RIEMANN_SOURCE)/rpt3_euler_mapgrid.f90 $(RIEMANN_SOURCE)/rptt3_euler_mapgrid.f90 +euler_3d_gmap.so: rpn3_euler_mapgrid.f90 rpt3_euler_mapgrid.f90 rptt3_euler_mapgrid.f90 f2py -m $(basename $(notdir $@)) -c $^ -mappedGrid.so: $(RIEMANN_SOURCE)/euler3d_mappedGrid.f90 +mappedGrid.so: euler3d_mappedGrid.f90 f2py -m $(basename $(notdir $@)) -c $^ clean: diff --git a/examples/euler_gravity_3d/__init__.py b/examples/euler_gravity_3d/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/examples/euler_gravity_3d/euler3d_mappedGrid.f90 b/examples/euler_gravity_3d/euler3d_mappedGrid.f90 new file mode 100644 index 000000000..55d149fdf --- /dev/null +++ b/examples/euler_gravity_3d/euler3d_mappedGrid.f90 @@ -0,0 +1,679 @@ +MODULE euler3d_mappedgrid + + REAL(kind=8), PARAMETER :: pi=3.1415926535897932385,twopi=2*pi + REAL(kind=8), PARAMETER :: rEarth=637120000.0 + +CONTAINS + !------------------------------------------------------------------------- + ! Name: + ! volumeOfHex(x,y,z): + ! + ! Description: + ! Compute Volume of an LD Hexahedron given the 8 vertices + ! from paper: + ! @book{Grandy_1997, + ! title={Efficient computation of volume of hexahedral cells}, + ! url={http://www.osti.gov/scitech/servlets/purl/632793}, + ! DOI={10.2172/632793}, + ! abstractNote={This report describes an efficient method to compute + ! the volume of hexahedral cells used in three-dimensional + ! hydrodynamics simulation. Two common methods for + ! creating the hexahedron using triangular boundaries are + ! considered.}, + ! author={Grandy, J.}, year={1997}, month={Oct}} + ! https://doi.org/10.2172/632793 + ! + ! 6*V_LD = + ! |x6-x0 x3-x0 x2-x7| |x6-x0 x4-x0 x7-x5| |x6-x0 x1-x0 x5-x2| + ! |y6-y0 y3-y0 y2-y7| + |y6-y0 y4-y0 y7-y5| + |y6-y0 y1-y0 y5-y2| + ! |z6-z0 z3-z0 z2-z7| |z6-z0 z4-z0 z7-z5| |z6-z0 z1-z0 z5-z2| + ! + ! Inputs: + ! x[:],y[:],z[:] : x,y, and z coordinates of eight hexahedron vertices + ! + ! Output: + ! volumeOfHex : Volume of LD hexahedron + !------------------------------------------------------------------------- + REAL(kind=8) FUNCTION volumeOfHex(x,y,z) + + IMPLICIT NONE + + REAL(kind=8), DIMENSION(8), INTENT(IN) :: x,y,z + + ! Compute Volume of Hexahedron + volumeOfHex = & + (x(7)-x(1))*((y(4)-y(1))*(z(3)-z(8)) - (y(3)-y(8))*(z(4)-z(1)))+& + (x(4)-x(1))*((y(3)-y(8))*(z(7)-z(1)) - (y(7)-y(1))*(z(3)-z(8)))+& + (x(3)-x(8))*((y(7)-y(1))*(z(4)-z(1)) - (y(4)-y(1))*(z(7)-z(1)))& + +& + (x(7)-x(1))*((y(5)-y(1))*(z(8)-z(6)) - (y(8)-y(6))*(z(5)-z(1)))+& + (x(5)-x(1))*((y(8)-y(6))*(z(7)-z(1)) - (y(7)-y(1))*(z(8)-z(6)))+& + (x(8)-x(6))*((y(7)-y(1))*(z(5)-z(1)) - (y(5)-y(1))*(z(7)-z(1)))& + +& + (x(7)-x(1))*((y(2)-y(1))*(z(6)-z(3)) - (y(6)-y(3))*(z(2)-z(1)))+& + (x(2)-x(1))*((y(6)-y(3))*(z(7)-z(1)) - (y(7)-y(1))*(z(6)-z(3)))+& + (x(6)-x(3))*((y(7)-y(1))*(z(2)-z(1)) - (y(2)-y(1))*(z(7)-z(1))) + volumeOfHex = ABS(volumeOfHex)/6.d0 + RETURN + END FUNCTION volumeOfHex + + !--------------------------------------------------------------------------- + ! Name: mapc2pWrapper(xc,yc,zc,mx,my,mz,xyz0,xyzN,mapType,xp,yp,zp) + ! + ! Description: mapping from computational to physical space + ! + ! Inputs: + ! xc,yc,zc : computational coordinate values + ! mx,my,mz : grid sizes + ! xyz0,xyzN : altitude at bottom and top of grid + ! mapType : type of mapping to compute + ! + ! Outputs: + ! xp,yp,zp : physical coordinate values + !--------------------------------------------------------------------------- + SUBROUTINE mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc + REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN + INTEGER, INTENT(IN) :: mz + CHARACTER(20) :: mapType + + ! Output + REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp + REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp + REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp + + ! Compute Physical Grid From Computational + SELECT CASE(mapType) + CASE("Identity") + CALL mapc2pIdentity(xc,yc,zc,xp,yp,zp) + CASE("ZeroToOne") + CALL mapc2pZeroToOne(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + CASE("Rotate") + CALL mapc2pRotate(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + CASE("Spherical") + CALL mapc2pSpherical(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + CASE DEFAULT + WRITE(*,*)"Invalid mapType: ",mapType + STOP + END SELECT + + RETURN + END SUBROUTINE mapc2pWrapper + + !--------------------------------------------------------------------------- + ! Name: mapc2pIdentity(xc,yc,zc,xp,yp,zp) + ! + ! Description: Identity mapping from computational to physical space + ! + ! Inputs: + ! xc,yc,zc : computational coordinate values + ! Outputs: + ! xp,yp,zp : physical coordinate values + !--------------------------------------------------------------------------- + SUBROUTINE mapc2pIdentity(xc,yc,zc,xp,yp,zp) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc + + ! Output + REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp + REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp + REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp + + ! x,xi coordinate + xp = xc + + ! y,eta coordinate + yp = yc + + ! z,phi coordinate + zp = zc + + RETURN + END SUBROUTINE mapc2pIdentity + + !--------------------------------------------------------------------------- + ! Name: mapc2pZeroToOne(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + ! + ! Description: mapping from computational to physical space using a + ! computational grid from 0 to 1 + ! + ! Inputs: + ! xc,yc,zc : computational coordinate values + ! xyz0,xyzN : physical coordinates at bottom and top of grid + ! + ! Outputs: + ! xp,yp,zp : physical coordinate values + !--------------------------------------------------------------------------- + SUBROUTINE mapc2pZeroToOne(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc + REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN + + ! Output + REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp + REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp + REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp + + ! x,xi coordinate + xp = xyz0(1) + xc*(xyzN(1)-xyz0(1)) + + ! y,eta coordinate + yp = xyz0(2) + yc*(xyzN(2)-xyz0(2)) + + ! z,phi coordinate + zp = xyz0(3) + zc*(xyzN(3)-xyz0(3)) + + RETURN + END SUBROUTINE mapc2pZeroToOne + + !--------------------------------------------------------------------------- + ! Name: mapc2pRotate(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + ! + ! Description: mapping from computational to physical space using a + ! computational grid from 0 to 1 and rotated by angle theta + ! + ! Inputs: + ! xc,yc,zc : computational coordinate values + ! xyz0,xyzN : physical coordinates at bottom and top of grid + ! + ! Outputs: + ! xp,yp,zp : physical coordinate values + !--------------------------------------------------------------------------- + SUBROUTINE mapc2pRotate(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc + REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN + + ! Output + REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp + REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp + REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp + + ! Local + REAL(kind=8) :: theta + REAL(kind=8), DIMENSION(SIZE(xc)) :: xp0 + REAL(kind=8), DIMENSION(SIZE(yc)) :: yp0 + REAL(kind=8), DIMENSION(SIZE(zc)) :: zp0 + + theta=0.d0*pi/180.d0 + + ! x,xi coordinate + xp = xyz0(1) + xc*(xyzN(1)-xyz0(1)) + + ! y,eta coordinate + yp = xyz0(2) + yc*(xyzN(2)-xyz0(2)) + + ! z,phi coordinate + zp = xyz0(3) + zc*(xyzN(3)-xyz0(3)) + + xp0 = xp + yp0 = yp + zp0 = zp + ! Rotate about x-axis by angle theta + xp = xp0 + yp = yp0*COS(theta) - zp0*SIN(theta) + zp = yp0*SIN(theta) + zp0*COS(theta) + + xp0 = xp + yp0 = yp + zp0 = zp + ! Rotate about y-axis by angle theta + xp = xp0*COS(theta) + zp0*SIN(theta) + yp = yp0 + zp =-xp0*SIN(theta) + zp0*COS(theta) + + xp0 = xp + yp0 = yp + zp0 = zp + ! Rotate about z-axis by angle theta + xp = xp0*COS(theta) - yp0*SIN(theta) + yp = xp0*SIN(theta) + yp0*COS(theta) + zp = zp0 + + RETURN + END SUBROUTINE mapc2pRotate + + !--------------------------------------------------------------------------- + ! Name: mapc2pSpherical(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + ! + ! Description: mapping from computational to physical space using a + ! computational grid from 0 to 1 and then using a + ! spherical transformation + ! + ! Inputs: + ! xc,yc,zc : computational coordinate values + ! xyz0,xyzN : physical coordinates at bottom and top of grid + ! + ! Outputs: + ! xp,yp,zp : physical coordinate values + !--------------------------------------------------------------------------- + SUBROUTINE mapc2pSpherical(xc,yc,zc,xyz0,xyzN,xp,yp,zp) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc + REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN + + ! Output + REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp + REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp + REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp + + ! Local + REAL(kind=8), DIMENSION(SIZE(xc)) :: theta + REAL(kind=8), DIMENSION(SIZE(yc)) :: phi + REAL(kind=8), DIMENSION(SIZE(zc)) :: r + + ! x,xi coordinate + xp = xyz0(1) + xc*(xyzN(1)-xyz0(1)) + + ! y,eta coordinate + yp = xyz0(2) + yc*(xyzN(2)-xyz0(2)) + + ! z,phi coordinate + zp = xyz0(3) + zc*(xyzN(3)-xyz0(3)) + + ! Spherical to Cartesian Coordinates + theta = xp + phi = yp + r = zp + xp = r*SIN(theta)*COS(phi) + yp = r*SIN(theta)*SIN(phi) + zp = r*COS(theta) + + RETURN + END SUBROUTINE mapc2pSpherical + + !--------------------------------------------------------------------------- + ! Name: determinant(A,detA) + ! + ! Description: Compute the determinant of a 3x3 array + ! + ! Inputs: A : 3x3 Array + ! + ! Outputs: detA : determinant of A + !--------------------------------------------------------------------------- + SUBROUTINE determinant(A,detA) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(3,3), INTENT(IN) :: A + + ! Output + REAL(kind=8), INTENT(OUT) :: detA + + detA = A(1,1)*A(2,2)*A(3,3) & + + A(2,1)*A(3,2)*A(1,3) & + + A(3,1)*A(1,2)*A(2,3) & + - A(3,1)*A(2,2)*A(1,3) & + - A(2,1)*A(1,2)*A(3,3) & + - A(1,1)*A(3,2)*A(2,3) + + RETURN + END SUBROUTINE determinant + + !--------------------------------------------------------------------------- + ! Name: unitNormal(a,b,c,n) + ! + ! Description: Compute the unit normal of given 3 indices in 3D space + ! + ! Inputs: a,b,c : indices with x,y, and z coordinates + ! + ! Outputs: n : unit normals of plane defined from points a,b,c + !--------------------------------------------------------------------------- + SUBROUTINE unitNormal(a,b,c,n) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(3), INTENT(IN) :: a,b,c + + ! Output + REAL(kind=8), DIMENSION(3), INTENT(OUT) :: n + + ! Local + REAL(kind=8), DIMENSION(3,3) :: D + REAL(kind=8) :: magnitude + + D(:,1) = (/1.d0,a(2),a(3)/) + D(:,2) = (/1.d0,b(2),b(3)/) + D(:,3) = (/1.d0,c(2),c(3)/) + CALL determinant(D,n(1)) + D(:,1) = (/a(1),1.d0,a(3)/) + D(:,2) = (/b(1),1.d0,b(3)/) + D(:,3) = (/c(1),1.d0,c(3)/) + CALL determinant(D,n(2)) + D(:,1) = (/a(1),a(2),1.d0/) + D(:,2) = (/b(1),b(2),1.d0/) + D(:,3) = (/c(1),c(2),1.d0/) + CALL determinant(D,n(3)) + + magnitude = SQRT(n(1)**2 + n(2)**2 + n(3)**2) + n(1) = n(1)/magnitude + n(2) = n(2)/magnitude + n(3) = n(3)/magnitude + + RETURN + END SUBROUTINE unitNormal + + !--------------------------------------------------------------------------- + ! Name: dot(a,b,adotb) + ! + ! Description: Compute the dot product of vectors a and b + ! + ! Inputs: a,b : 3D vectors + ! + ! Outputs: adotb : dot product of a and b + !--------------------------------------------------------------------------- + SUBROUTINE dot(a,b,adotb) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(3), INTENT(IN) :: a,b + + ! Output + REAL(kind=8), INTENT(OUT) :: adotb + + adotb = a(1)*b(1) + a(2)*b(2) + a(3)*b(3) + + RETURN + END SUBROUTINE dot + + !--------------------------------------------------------------------------- + ! Name: cross(a,b,c) + ! + ! Description: Compute the cross product of vectors a and b + ! + ! Inputs: a,b : 3D vectors + ! + ! Outputs: c : cross product of a and b + !--------------------------------------------------------------------------- + SUBROUTINE cross(a,b,c) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(3), INTENT(IN) :: a,b + + ! Output + REAL(kind=8), DIMENSION(3), INTENT(OUT) :: c + + c(1) = a(2)*b(3) - a(3)*b(2) + c(2) = a(3)*b(1) - a(1)*b(3) + c(3) = a(1)*b(2) - a(2)*b(1) + + RETURN + END SUBROUTINE cross + + !--------------------------------------------------------------------------- + ! Name: areaPolygon(poly) + ! + ! Description: Compute the area of a polygon given its vertices + ! + ! Inputs: poly : vertices of polygon + ! + ! Outputs: area : area of polygon + !--------------------------------------------------------------------------- + SUBROUTINE areaPolygon(poly,polyArea) + + IMPLICIT NONE + + ! Input + REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: poly + + ! Output + REAL(kind=8), INTENT(OUT) :: polyArea + + ! Local + INTEGER :: i,np + REAL(kind=8), DIMENSION(3) :: total,prod,vi1,vi2,n + + ! Compute Area of a Polygon Given its Vertices + np = SIZE(poly,2) + IF(np < 3)THEN + WRITE(*,*)"This is not a polygon, no area" + STOP + END IF + total=0 + DO i = 1,np + vi1 = poly(:,i) + IF(i == np)THEN + vi2 = poly(:,1) + ELSE + vi2 = poly(:,i+1) + END IF + CALL cross(vi1,vi2,prod) + total(1) = total(1) + prod(1) + total(2) = total(2) + prod(2) + total(3) = total(3) + prod(3) + END DO + CALL unitNormal(poly(:,1),poly(:,2),poly(:,3),n) + CALL dot(total,n,polyArea) + polyArea = ABS(polyArea*0.5d0) + + RETURN + END SUBROUTINE areaPolygon + + !--------------------------------------------------------------------------- + ! Name: setAuxiliaryVariables(num_aux,mbc,mx,my,mz,xlower,ylower,zlower,& + ! dxc,dyc,dzc,xyz0,xyzN,aux) + ! + ! Description: Compute auxiliary variables for non-uniform/mapped grids + ! + ! Inputs: + ! xc,yc,zc : computational coordinate values + ! mx,my,mz : grid sizes + ! xyz0,xyzN : altitude at bottom and top of grid + ! + ! Outputs: + ! xp,yp,zp : physical coordinate values + !--------------------------------------------------------------------------- + SUBROUTINE setAuxiliaryVariables(num_aux,mbc,mx,my,mz,xlower,ylower,zlower,& + dxc,dyc,dzc,xyz0,xyzN,mapType,aux) + + IMPLICIT NONE + + ! Input + INTEGER, INTENT(IN) :: num_aux,mbc,mx,my,mz + REAL(kind=8), INTENT(IN) :: xlower,ylower,zlower,dxc,dyc,dzc + REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN + CHARACTER(20) :: mapType + + ! Output + REAL(kind=8), DIMENSION(num_aux,mx+2*mbc,my+2*mbc,mz+2*mbc), & + INTENT(OUT) :: aux + + ! Local + REAL(kind=8), DIMENSION(3) :: norm + REAL(kind=8), DIMENSION(3,4) :: poly + REAL(kind=8), DIMENSION(8) :: xcCorner,ycCorner,zcCorner,& + xpCorner,ypCorner,zpCorner + INTEGER :: i,j,k + REAL(kind=8) :: vHex,area + REAL(kind=8), DIMENSION(1) :: xc,yc,zc,xp,yp,zp + + ! Compute Auxiliary Variables + DO i = 1,mx+2*mbc + DO j = 1,my+2*mbc + DO k = 1,mz+2*mbc + + ! Compute Physical Coordinates of Corners of Hexahedron + + ! i-1/2, j-1/2, k-1/2 Corner + xcCorner(1) = xlower + (i-1)*dxc + ycCorner(1) = ylower + (j-1)*dyc + zcCorner(1) = zlower + (k-1)*dzc + xc(1) = xcCorner(1) + yc(1) = ycCorner(1) + zc(1) = zcCorner(1) + CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + xpCorner(1) = xp(1) + ypCorner(1) = yp(1) + zpCorner(1) = zp(1) + + ! i+1/2, j-1/2, k-1/2 Corner + xcCorner(2) = xcCorner(1) + dxc + ycCorner(2) = ycCorner(1) + zcCorner(2) = zcCorner(1) + xc(1) = xcCorner(2) + yc(1) = ycCorner(2) + zc(1) = zcCorner(2) + CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + xpCorner(2) = xp(1) + ypCorner(2) = yp(1) + zpCorner(2) = zp(1) + + ! i+1/2, j+1/2, k-1/2 Corner + xcCorner(3) = xcCorner(1) + dxc + ycCorner(3) = ycCorner(1) + dyc + zcCorner(3) = zcCorner(1) + xc(1) = xcCorner(3) + yc(1) = ycCorner(3) + zc(1) = zcCorner(3) + CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + xpCorner(3) = xp(1) + ypCorner(3) = yp(1) + zpCorner(3) = zp(1) + + ! i-1/2, j+1/2, k-1/2 Corner + xcCorner(4) = xcCorner(1) + ycCorner(4) = ycCorner(1) + dyc + zcCorner(4) = zcCorner(1) + xc(1) = xcCorner(4) + yc(1) = ycCorner(4) + zc(1) = zcCorner(4) + CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + xpCorner(4) = xp(1) + ypCorner(4) = yp(1) + zpCorner(4) = zp(1) + + ! i-1/2, j-1/2, k+1/2 Corner + xcCorner(5) = xcCorner(1) + ycCorner(5) = ycCorner(1) + zcCorner(5) = zcCorner(1) + dzc + xc(1) = xcCorner(5) + yc(1) = ycCorner(5) + zc(1) = zcCorner(5) + CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + xpCorner(5) = xp(1) + ypCorner(5) = yp(1) + zpCorner(5) = zp(1) + + ! i+1/2, j-1/2, k+1/2 Corner + xcCorner(6) = xcCorner(1) + dxc + ycCorner(6) = ycCorner(1) + zcCorner(6) = zcCorner(1) + dzc + xc(1) = xcCorner(6) + yc(1) = ycCorner(6) + zc(1) = zcCorner(6) + CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + xpCorner(6) = xp(1) + ypCorner(6) = yp(1) + zpCorner(6) = zp(1) + + ! i+1/2, j+1/2, k+1/2 Corner + xcCorner(7) = xcCorner(1) + dxc + ycCorner(7) = ycCorner(1) + dyc + zcCorner(7) = zcCorner(1) + dzc + xc(1) = xcCorner(7) + yc(1) = ycCorner(7) + zc(1) = zcCorner(7) + CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + xpCorner(7) = xp(1) + ypCorner(7) = yp(1) + zpCorner(7) = zp(1) + + ! i-1/2, j+1/2, k+1/2 Corner + xcCorner(8) = xcCorner(1) + ycCorner(8) = ycCorner(1) + dyc + zcCorner(8) = zcCorner(1) + dzc + xc(1) = xcCorner(8) + yc(1) = ycCorner(8) + zc(1) = zcCorner(8) + CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) + xpCorner(8) = xp(1) + ypCorner(8) = yp(1) + zpCorner(8) = zp(1) + + ! Compute Aux Variables (Normals/Area Ratios/Volume Ratio) + + ! \xi normal + poly(:,1) = (/xpCorner(1),ypCorner(1),zpCorner(1)/) + poly(:,2) = (/xpCorner(4),ypCorner(4),zpCorner(4)/) + poly(:,3) = (/xpCorner(8),ypCorner(8),zpCorner(8)/) + poly(:,4) = (/xpCorner(5),ypCorner(5),zpCorner(5)/) + CALL unitNormal(poly(:,1),poly(:,2),poly(:,3),norm) + aux(1,i,j,k) = norm(1) ! nx(i-1/2,j,k) + aux(2,i,j,k) = norm(2) ! ny(i-1/2,j,k) + aux(3,i,j,k) = norm(3) ! nz(i-1/2,j,k) + CALL areaPolygon(poly,area) + aux(4,i,j,k) = area/(dyc*dzc) ! gamma(i-1/2,j,k) + + ! \eta normal + poly(:,1) = (/xpCorner(1),ypCorner(1),zpCorner(1)/) + poly(:,2) = (/xpCorner(5),ypCorner(5),zpCorner(5)/) + poly(:,3) = (/xpCorner(6),ypCorner(6),zpCorner(6)/) + poly(:,4) = (/xpCorner(2),ypCorner(2),zpCorner(2)/) + CALL unitNormal(poly(:,1),poly(:,2),poly(:,3),norm) + aux(5,i,j,k) = norm(1) ! nx(i,j-1/2,k) + aux(6,i,j,k) = norm(2) ! ny(i,j-1/2,k) + aux(7,i,j,k) = norm(3) ! nz(i,j-1/2,k) + CALL areaPolygon(poly,area) + aux(8,i,j,k) = area/(dxc*dzc) ! gamma(i,j-1/2,k) + + ! \phi normal + poly(:,1) = (/xpCorner(1),ypCorner(1),zpCorner(1)/) + poly(:,2) = (/xpCorner(2),ypCorner(2),zpCorner(2)/) + poly(:,3) = (/xpCorner(3),ypCorner(3),zpCorner(3)/) + poly(:,4) = (/xpCorner(4),ypCorner(4),zpCorner(4)/) + CALL unitNormal(poly(:,1),poly(:,2),poly(:,3),norm) + aux(9,i,j,k) = norm(1) ! nx(i,j,k-1/2) + aux(10,i,j,k) = norm(2) ! ny(i,j,k-1/2) + aux(11,i,j,k) = norm(3) ! nz(i,j,k-1/2) + CALL areaPolygon(poly,area) + aux(12,i,j,k) = area/(dxc*dyc) ! gamma(i,j,k-1/2) + + ! Volume of hexahedron + vHex = volumeOfHex(xpCorner,ypCorner,zpCorner) + + ! Capacity Function (Ratio of Physical:Computational volume) + aux(13,i,j,k) = vHex/(dxc*dyc*dzc) ! kappa(i,j,k) + + ! Grid Delta in Direction of Gravity + SELECT CASE(mapType) + CASE("Identity","ZeroToOne","Rotate") + aux(14,i,j,k) = ABS(zpCorner(5)-zpCorner(1)) ! dz + aux(15,i,j,k) = (zpCorner(5) + zpCorner(1))/2.d0 ! z + CASE("Spherical") + aux(14,i,j,k) = & + ABS(SQRT(xpCorner(5)**2+ypCorner(5)**2+zpCorner(5)**2) & + - SQRT(xpCorner(1)**2+ypCorner(1)**2+zpCorner(1)**2)) ! dr + aux(15,i,j,k) = (SQRT(xpCorner(5)**2+ypCorner(5)**2+zpCorner(5)**2) & + + SQRT(xpCorner(1)**2+ypCorner(1)**2+zpCorner(1)**2))/2.d0 & + - rEarth + END SELECT + + END DO + END DO + END DO + + RETURN + END SUBROUTINE setAuxiliaryVariables + +END MODULE euler3d_mappedgrid diff --git a/examples/euler_gravity_3d/rising_hot_sphere.py b/examples/euler_gravity_3d/rising_hot_sphere.py index d8d07923f..24e54a330 100755 --- a/examples/euler_gravity_3d/rising_hot_sphere.py +++ b/examples/euler_gravity_3d/rising_hot_sphere.py @@ -11,13 +11,15 @@ density (rho), x,y, and z momentum (rho*u,rho*v,rho*w), and energy. """ import numpy as np -from mappedGrid import euler3d_mappedgrid as mg +from clawpack.pyclaw.examples.euler_gravity_3d.mappedGrid import euler3d_mappedgrid as mg +# from .mappedGrid import euler3d_mappedgrid as mg try: from mpi4py import MPI mpiAvailable = True except ImportError: - raise ImportError('mpi4py is not available') + import warnings + warnings.warn('mpi4py is not available') mpiAvailable = False if mpiAvailable: @@ -175,8 +177,8 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ import logging solver.logger.setLevel(logging.DEBUG) - - import euler_3d_gmap + from clawpack.pyclaw.examples.euler_gravity_3d import euler_3d_gmap + # from . import euler_3d_gmap solver.rp = euler_3d_gmap solver.num_eqn = 5 solver.num_waves = 3 @@ -340,7 +342,7 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ claw.solver = solver claw.output_format = output_format claw.output_file_prefix = file_prefix - claw.keep_copy = False + claw.keep_copy = True if disable_output: claw.output_format = None claw.tfinal = tfinal diff --git a/examples/euler_gravity_3d/rising_hot_sphere_spherical.py b/examples/euler_gravity_3d/rising_hot_sphere_spherical.py index 88e7ac3a7..41ce20c08 100755 --- a/examples/euler_gravity_3d/rising_hot_sphere_spherical.py +++ b/examples/euler_gravity_3d/rising_hot_sphere_spherical.py @@ -11,14 +11,16 @@ density (rho), x,y, and z momentum (rho*u,rho*v,rho*w), and energy. """ import numpy as np -from mappedGrid import euler3d_mappedgrid as mg +from clawpack.pyclaw.examples.euler_gravity_3d.mappedGrid import euler3d_mappedgrid as mg +# from .mappedGrid import euler3d_mappedgrid as mg # Test for MPI, and set sizes accordingly try: from mpi4py import MPI mpiAvailable = True except ImportError: - raise ImportError('mpi4py is not available') + import warnings + warnings.warn('mpi4py is not available') mpiAvailable = False if mpiAvailable: @@ -172,8 +174,8 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ import logging solver.logger.setLevel(logging.DEBUG) - - import euler_3d_gmap + from clawpack.pyclaw.examples.euler_gravity_3d import euler_3d_gmap + #from . import euler_3d_gmap solver.rp = euler_3d_gmap solver.num_eqn = 5 solver.num_waves = 3 diff --git a/examples/euler_gravity_3d/rpn3_euler_mapgrid.f90 b/examples/euler_gravity_3d/rpn3_euler_mapgrid.f90 new file mode 100644 index 000000000..026b16998 --- /dev/null +++ b/examples/euler_gravity_3d/rpn3_euler_mapgrid.f90 @@ -0,0 +1,312 @@ +!----------------------------------------------------------------------------------- +! Name: rpn3(ixyz,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq) +! +! Description: Roe f-wave solver for the 3D Euler equations in general geometries +! +! Waves: 3 +! Equations: 5 +! +! Conserved Quantities: +! 1 density +! 2 x_momentum +! 3 y_momentum +! 4 z_momentum +! 5 energy +! +! Auxilliary Quantities: +! 1 nx xface (i-1/2,j,k) +! 2 ny xface (i-1/2,j,k) +! 3 nz xface (i-1/2,j,k) +! 4 gamma_x (area/(dyc*dzc)) +! 5 nx yface (i,j-1/2,k) +! 6 ny yface (i,j-1/2,k) +! 7 nz yface (i,j-1/2,k) +! 8 gamma_x (area/(dxc*dzc)) +! 9 nx zface (i,j,k-1/2) +! 10 ny zface (i,j,k-1/2) +! 11 nz zface (i,j,k-1/2) +! 12 gamma_z (area/(dxc*dyc)) +! 13 kappa(i,j,k) volumeOfHex/(dxc*dyc*dzc) +! 14 dz (Cartesian), dr (Spherical) +! 15 z (Cartesian), r (Spherical) +! +! Inputs: +! ixyz : direction to take slice x-direction if ixyz=1 +! y-direction if ixyz=2. +! z-direction if ixyz=3. +! maxm : max number of grid cells (less the ghost cells) +! meqn : number of equations in the system (=5) +! mwaves : nuber of waves in the system (=3) +! maux : number of auxilary variables (=15) +! mbc : number of ghost cells on either end +! mx : number of elements +! ql : state vector at left edge of each cell +! Note the i'th Riemann problem has left state qr(:,i-1) +! qr : state vector at right edge of each cell +! Note the i'th Riemann problem has right state ql(:,i) +! auxl : state of auxilary variable on left egde of cell +! auxr : state of auxilary variable on right egde of cell +! +! Outputs: +! fwave : f-wave vectors (eigenvectors of Roe matrix) +! s : eigenvalues (eigenvalues of Roe matrix) +! amdq : left-going fluctuations +! apdq : right-going fluctuations +! +! Adapted from rpn3_euler.f90 in $CLAWHOME/riemann/src +!----------------------------------------------------------------------------------- +SUBROUTINE rpn3(ixyz,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq) + + IMPLICIT NONE + + ! Input + INTEGER, INTENT(IN) :: ixyz,maxm,meqn,mwaves,maux,mbc,mx + REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(IN) :: ql,qr + REAL(kind=8), DIMENSION(maux,1-mbc:maxm+mbc), INTENT(IN) :: auxl,auxr + + ! Output + REAL(kind=8), DIMENSION(meqn,mwaves,1-mbc:maxm+mbc), INTENT(INOUT) :: fwave + REAL(kind=8), DIMENSION(mwaves,1-mbc:maxm+mbc), INTENT(INOUT) :: s + REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(INOUT) :: amdq,apdq + + ! Local Storage + REAL(kind=8), DIMENSION(meqn) :: delta,beta + REAL(kind=8), DIMENSION(-1:maxm+mbc) :: u,v,w,enth,ek,a,a2g1,psi + REAL(kind=8) :: gamma1,asqrd,pl,pr,vnl,vnr,rhsq2,rhsqrtl,rhsqrtr,nx,ny,nz,vn,& + areaRatio,dr,finterp,qinterp,quinterp,qvinterp,qwinterp + INTEGER :: i,j,mu=2,mv=3,mw=4,mws,nu,nv,nw,iratio=4,gFlux=0 + CHARACTER(20) :: nonsingularNormal + + ! Common block storage for ideal gas constant and gravity term + REAL(kind=8) :: gamma ! Ideal gas constant + REAL(kind=8) :: g_r ! Gravitational Constant Magnitude + LOGICAL :: gravity,gravityflux ! Turn Gravity Source Term On/Off + COMMON /cparam/ gamma,g_r,gravity,gravityflux + + ! Set (gamma-1) + gamma1 = gamma - 1.d0 + + ! Gravity in Energy Flux On/Off + gFlux = 0 + IF(gravityflux)gFlux=1 + + ! Make sure the number of waves is 3 for this solver + IF (mwaves /= 3) THEN + WRITE(6,*) '*** Must set mwaves=3 for this Riemann solver' + STOP + ENDIF + + ! Set mu to point to the component of the system that corresponds to momentum in + ! the direction of this slice, mv and mw to the orthogonal momentums. Set the + ! nu,nv,and nw indices for the normal vectors in a similar way. The iratio is + ! the aux array index for the gamma area scaling for mapped grids. + IF(ixyz == 1)THEN + mu = 2 + mv = 3 + mw = 4 + iratio = 4 + ELSE IF(ixyz == 2)THEN + mu = 3 + mv = 4 + mw = 2 + iratio = 8 + ELSE IF (ixyz == 3)THEN + mu = 4 + mv = 2 + mw = 3 + iratio = 12 + ELSE + WRITE(*,*)"Something is wrong with ixyz..." + STOP + ENDIF + nu = (ixyz-1)*4 + mu - 1 + nv = (ixyz-1)*4 + mv - 1 + nw = (ixyz-1)*4 + mw - 1 + + ! Note that notation for u,v, and w reflects assumption that the + ! Riemann problems are in the x-direction with u in the normal + ! direction and v and w in the orthogonal directions, but with the + ! above definitions of mu, mv, and mw the routine also works with + ! ixyz=2 and ixyz = 3 + ! and returns, for example, f0 as the Godunov flux g0 for the + ! Riemann problems u_t + g(u)_y = 0 in the y-direction. + + ! Initialize fwaves and speeds + fwave = 0.d0 + s = 0.d0 + + ! Loop over grid cell interfaces + DO i = 2-mbc, mx+mbc + IF (qr(1,i-1) <= 0.d0 .OR. ql(1,i) <= 0.d0) THEN + WRITE(*,*) i, mu, mv, mw + WRITE(*,'(5e12.4)') (qr(j,i-1),j=1,5) + WRITE(*,'(5e12.4)') (ql(j,i),j=1,5) + IF (ixyz == 1) WRITE(6,*) '*** (rpn3) rho <= 0 in x-sweep at ',i + IF (ixyz == 2) WRITE(6,*) '*** (rpn3) rho <= 0 in y-sweep at ',i + IF (ixyz == 3) WRITE(6,*) '*** (rpn3) rho <= 0 in z-sweep at ',i + WRITE(6,*) 'stopped with rho <= 0...' + STOP + ENDIF + rhsqrtl = SQRT(qr(1,i-1)) + rhsqrtr = SQRT(ql(1,i)) + pl = gamma1*(qr(5,i-1) - 0.5d0*(qr(mu,i-1)**2 + & + qr(mv,i-1)**2 + qr(mw,i-1)**2)/qr(1,i-1) & + - qr(1,i-1)*(g_r)*auxr(15,i-1)*gFlux ) + pr = gamma1*(ql(5,i) - 0.5d0*(ql(mu,i)**2 + & + ql(mv,i)**2 + ql(mw,i)**2)/ql(1,i) & + - ql(1,i)*(g_r)*auxl(15,i)*gFlux ) + rhsq2 = rhsqrtl + rhsqrtr + u(i) = (qr(mu,i-1)/rhsqrtl + ql(mu,i)/rhsqrtr) / rhsq2 + v(i) = (qr(mv,i-1)/rhsqrtl + ql(mv,i)/rhsqrtr) / rhsq2 + w(i) = (qr(mw,i-1)/rhsqrtl + ql(mw,i)/rhsqrtr) / rhsq2 + enth(i) = ((qr(5,i-1)+pl)/rhsqrtl + (ql(5,i)+pr)/rhsqrtr) / rhsq2 + psi(i) = (qr(1,i-1)*(g_r)*auxr(15,i-1)/rhsqrtl & + + ql(1,i)*(g_r)*auxl(15,i)/rhsqrtr)/rhsq2*gFlux + ek(i) = 0.5d0*(u(i)**2 + v(i)**2 + w(i)**2) + asqrd = gamma1*(enth(i) - ek(i)) + IF (i>=0 .AND. i<=mx .AND. asqrd <= 0.d0) THEN + IF (ixyz == 1) WRITE(6,*) '*** (rpn3) a**2 <= 0 in x-sweep at ',i + IF (ixyz == 2) WRITE(6,*) '*** (rpn3) a**2 <= 0 in y-sweep at ',i + IF (ixyz == 3) WRITE(6,*) '*** (rpn3) a**2 <= 0 in z-sweep at ',i + WRITE(6,*) 'stopped with a**2 < 0...' + STOP + ENDIF + a(i) = SQRT(asqrd) + a2g1(i) = asqrd/gamma1 + END DO + + ! Now split the jump in q1d at each interface into f-waves + ! find alpha(1) thru alpha(5), the coefficients of the 5 eigenvectors: + DO i = 2-mbc, mx+mbc + + ! Normal Vectors, Normal Velocity, and Area Ratio + nx = auxl(nu,i) + ny = auxl(nv,i) + nz = auxl(nw,i) + vn = u(i)*nx + v(i)*ny + w(i)*nz + areaRatio = auxl(iratio,i) + + ! Compute the jumps (delta F) in the flux functions + pl = gamma1*(ql(5,i) & + - 0.5d0*(ql(mu,i)**2+ql(mv,i)**2+ql(mw,i)**2)/ql(1,i) & + - ql(1,i)*(g_r)*auxl(15,i)*gFlux) + pr = gamma1*(qr(5,i-1) & + - 0.5d0*(qr(mu,i-1)**2+qr(mv,i-1)**2+qr(mw,i-1)**2)/qr(1,i-1) & + - qr(1,i-1)*(g_r)*auxr(15,i-1)*gFlux) + vnl = (ql(mu,i)*nx + ql(mv,i)*ny + ql(mw,i)*nz)/ql(1,i) + vnr = (qr(mu,i-1)*nx + qr(mv,i-1)*ny + qr(mw,i-1)*nz)/qr(1,i-1) + delta(1) = ql(1,i)*vnl - qr(1,i-1)*vnr + delta(2) = (vnl*ql(mu,i) + pl*nx) - (vnr*qr(mu,i-1) + pr*nx) + delta(3) = (vnl*ql(mv,i) + pl*ny) - (vnr*qr(mv,i-1) + pr*ny) + delta(4) = (vnl*ql(mw,i) + pl*nz) - (vnr*qr(mw,i-1) + pr*nz) + delta(5) = ((ql(5,i) + pl)*vnl) - ((qr(5,i-1) + pr)*vnr) + + ! Modify delta(2/3/4/5) to account for gravitational term + IF(gravity .AND. ixyz==3)THEN + dr = (auxl(14,i) + auxl(14,i-1))*0.5d0 + finterp = auxl(14,i-1)/(auxl(14,i)+auxl(14,i-1)) + qinterp = ql(1,i)*finterp + qr(1,i-1)*(1.d0-finterp) + quinterp = ql(mu,i)*finterp + qr(mu,i-1)*(1.d0-finterp) + qvinterp = ql(mv,i)*finterp + qr(mv,i-1)*(1.d0-finterp) + qwinterp = ql(mw,i)*finterp + qr(mw,i-1)*(1.d0-finterp) + delta(2) = delta(2) + dr*g_r*nx*qinterp + delta(3) = delta(3) + dr*g_r*ny*qinterp + delta(4) = delta(4) + dr*g_r*nz*qinterp + IF(.NOT.gravityflux)& + delta(5) = delta(5) + dr*g_r*(quinterp*nx + qvinterp*ny + qwinterp*nz) + END IF + + ! Pick the largest normal to avoid using a singular direction + IF(ABS(nx) >= ABS(ny) .AND. ABS(nx) >= ABS(nz))THEN + nonsingularNormal = "x" + ELSEIF(ABS(ny) >= ABS(nx) .AND. ABS(ny) >= ABS(nz))THEN + nonsingularNormal = "y" + ELSEIF(ABS(nz) >= ABS(nx) .AND. ABS(nz) >= ABS(ny))THEN + nonsingularNormal = "z" + ELSE + WRITE(*,*)"Invalid Normal Vector..." + WRITE(*,*)"nx: ",nx,"ny: ",ny,"nz: ",nz + STOP + END IF + + beta(4) = ((enth(i) + psi(i) - 2.d0*ek(i))*delta(1) & + + u(i)*delta(2) + v(i)*delta(3) + w(i)*delta(4) - delta(5))/a2g1(i) + beta(5) = ( (a(i)-vn)*delta(1) - a(i)*beta(4) & + + nx*delta(2) + ny*delta(3) + nz*delta(4) )/(2.d0*a(i)) + beta(1) = delta(1) - beta(4) - beta(5) + + SELECT CASE(TRIM(nonsingularNormal)) + CASE("x") + beta(2) = ( (v(i)-vn*ny)*delta(1) + nx*ny*delta(2) & + + (ny**2-1.d0)*delta(3) + ny*nz*delta(4) )/nx + beta(3) = (-(w(i)-vn*nz)*delta(1) - nx*nz*delta(2) & + - ny*nz*delta(3) - (nz**2-1.d0)*delta(4) )/nx + CASE("y") + beta(2) = (-(u(i)-vn*nx)*delta(1) - (nx**2-1.d0)*delta(2) & + - ny*nx*delta(3) - nx*nz*delta(4) )/ny + beta(3) = ( (w(i)-vn*nz)*delta(1) + nx*nz*delta(2) & + + ny*nz*delta(3) + (nz**2-1.d0)*delta(4) )/ny + CASE("z") + beta(2) = ( (u(i)-vn*nx)*delta(1) + (nx**2-1.d0)*delta(2) & + + nx*ny*delta(3) + nz*nx*delta(4) )/nz + beta(3) = (-(v(i)-vn*ny)*delta(1) - nx*ny*delta(2) & + - (ny**2-1.d0)*delta(3) - nz*ny*delta(4) )/nz + END SELECT + + ! Compute the f-waves. + ! Note that the 2-wave, 3-wave and 4-wave travel at the same speed + ! and are lumped together in fwave(.,2,.). The 5-wave is then stored + ! in fwave(.,3,.). + fwave(1,1,i) = beta(1) + fwave(mu,1,i) = beta(1)*(u(i)-a(i)*nx) + fwave(mv,1,i) = beta(1)*(v(i)-a(i)*ny) + fwave(mw,1,i) = beta(1)*(w(i)-a(i)*nz) + fwave(5,1,i) = beta(1)*(enth(i) + psi(i) - a(i)*vn) + s(1,i) = (vn-a(i))*areaRatio + + fwave(1,2,i) = beta(4) + SELECT CASE(TRIM(nonsingularNormal)) + CASE("x") + fwave(mu,2,i) = beta(4)*u(i) + beta(2)*ny - beta(3)*nz + fwave(mv,2,i) = beta(4)*v(i) - beta(2)*nx + fwave(mw,2,i) = beta(4)*w(i) + beta(3)*nx + fwave(5,2,i) = beta(4)*(ek(i) + psi(i)) & + + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(w(i)*nx-u(i)*nz) + CASE("y") + fwave(mu,2,i) = beta(4)*u(i) + beta(2)*ny + fwave(mv,2,i) = beta(4)*v(i) - beta(2)*nx + beta(3)*nz + fwave(mw,2,i) = beta(4)*w(i) - beta(3)*ny + fwave(5,2,i) = beta(4)*(ek(i) + psi(i)) & + + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(v(i)*nz-w(i)*ny) + CASE("z") + fwave(mu,2,i) = beta(4)*u(i) - beta(2)*nz + fwave(mv,2,i) = beta(4)*v(i) + beta(3)*nz + fwave(mw,2,i) = beta(4)*w(i) + beta(2)*nx - beta(3)*ny + fwave(5,2,i) = beta(4)*(ek(i) + psi(i)) & + + beta(2)*(w(i)*nx-u(i)*nz) + beta(3)*(v(i)*nz-w(i)*ny) + END SELECT + s(2,i) = vn*areaRatio + + fwave(1,3,i) = beta(5) + fwave(mu,3,i) = beta(5)*(u(i)+a(i)*nx) + fwave(mv,3,i) = beta(5)*(v(i)+a(i)*ny) + fwave(mw,3,i) = beta(5)*(w(i)+a(i)*nz) + fwave(5,3,i) = beta(5)*(enth(i) + psi(i) + a(i)*vn) + s(3,i) = (vn+a(i))*areaRatio + + ! Compute fluctuations amdq and apdq + ! amdq = SUM s*fwave over left-going waves + ! apdq = SUM s*fwave over right-going waves + amdq(:,i) = 0.d0 + apdq(:,i) = 0.d0 + DO mws=1,mwaves + fwave(:,mws,i) = fwave(:,mws,i)*areaRatio + IF (s(mws,i) < 0.d0) THEN + amdq(:,i) = amdq(:,i) + fwave(:,mws,i) + ELSE + apdq(:,i) = apdq(:,i) + fwave(:,mws,i) + ENDIF + END DO + + END DO +END SUBROUTINE rpn3 diff --git a/examples/euler_gravity_3d/rpt3_euler_mapgrid.f90 b/examples/euler_gravity_3d/rpt3_euler_mapgrid.f90 new file mode 100644 index 000000000..09f215eaf --- /dev/null +++ b/examples/euler_gravity_3d/rpt3_euler_mapgrid.f90 @@ -0,0 +1,295 @@ +!----------------------------------------------------------------------------------- +! Name: rpt3(ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq, +! bmasdq,bpasdq) +! +! Description: Roe f-wave solver in the transverse direction for the +! 3D Euler equations in general geometries +! +! Inputs: +! ixyz : direction to take slice x-direction if ixyz=1 +! y-direction if ixyz=2. +! z-direction if ixyz=3. +! icoor : direction in which the transverse solve should be +! performed. icoor=2 (split in y-like direction) +! icoor=3 (split in z-like direction) +! ixyz=1, icoor=2 means asdq=A^*\Dq, split in y into: +! bmasdq = B^-A^*\Dq, +! bpasdq = B^+A^*\Dq. +! ixyz=2, icoor=2 means asdq=B^*\Dq, split in z into: +! bmasdq = C^-B^*\Dq, +! bpasdq = C^+B^*\Dq. +! imp : index for aux arrays +! maxm : max number of grid cells (less the ghost cells) +! meqn : number of equations in the system +! mwaves : nuber of waves in the system +! maux : number of auxilary equations +! mbc : number of ghost cells on either end +! mx : number of elements +! ql : state vector at left edge of each cell +! Note the i'th Riemann problem has left state qr(:,i-1) +! qr : state vector at right edge of each cell +! Note the i'th Riemann problem has right state ql(:,i) +! aux1 : +! aux2 : +! aux3 : +! asdq : array of flux differences (A*\Dq). asdq(:,i) is the +! flux difference propagating away from the interface +! between cells i-1 and i. Note: asdq represents +! B*\Dq if ixyz=2 or C*\Dq if ixyz=3. +! +! Outputs: +! bmasdq : left-going flux differences +! bpasdq : right-going flux differences +! +! Adapted from rpt3_euler.f90 in $CLAWHOME/riemann/src +!----------------------------------------------------------------------------------- +SUBROUTINE rpt3(ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,& + asdq,bmasdq,bpasdq) + + IMPLICIT NONE + + ! Input + INTEGER, INTENT(IN) :: ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx + REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(IN) :: ql,qr,asdq + REAL(kind=8), DIMENSION(maux,1-mbc:maxm+mbc,3), INTENT(IN) :: aux1,aux2,aux3 + + ! Output + REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(INOUT) :: bmasdq,bpasdq + + ! Local Storage + REAL(kind=8), DIMENSION(meqn) :: beta + REAL(kind=8), DIMENSION(-1:maxm+mbc) :: u2v2w2,u,v,w,enth,a,g1a2,euv + REAL(kind=8) :: gamma1,asqrd,pl,pr,rhsq2,rhsqrtl,rhsqrtr,nx,ny,nz,vn,& + areaRatioMinus,areaRatioPlus + REAL(kind=8), DIMENSION(meqn,mwaves) :: waveb + REAL(kind=8), DIMENSION(mwaves) :: sb + INTEGER :: i,j,mu=2,mv=3,mw=4,mws,nu=1,nv=2,nw=3,iratio=4 + CHARACTER(20) :: nonsingularNormal + + ! Common block storage for ideal gas constant + REAL(kind=8) :: gamma + COMMON /cparam/ gamma + + ! Set (gamma-1) + gamma1 = gamma - 1.d0 + + ! Make sure the number of waves is 3 for the 3D Euler equations + IF (mwaves /= 3) THEN + WRITE(6,*) '*** Should have mwaves=3 for this Riemann solver' + STOP + ENDIF + + ! Set mu to point to the component of the system that corresponds to momentum in + ! the direction of this slice, mv and mw to the orthogonal momentums. + ! + ! ixyz indicates the direction of the original Riemann solve, + ! called the x-like direction in the table below: + ! + ! x-like direction y-like direction z-like direction + ! ixyz=1: x y z + ! ixyz=2: y z x + ! ixyz=3: z x y + IF(ixyz == 1)THEN ! x-like direction + mu = 2 + mv = 3 + mw = 4 + ! Choose the coordinate direction to solve the Riemann problem + ! Permute nu,nv,nw for the original ixyz incident direction (x-like) + IF(icoor == 2)THEN ! y-like direction + nu = 5 + nv = 6 + nw = 7 + iratio = 8 + ELSEIF(icoor == 3)THEN ! z-like direction + nu = 9 + nv = 10 + nw = 11 + iratio = 12 + END IF + ELSEIF(ixyz == 2)THEN ! y-like direction + mu = 3 + mv = 4 + mw = 2 + ! Choose the coordinate direction to solve the Riemann problem + ! Permute nu,nv,nw for the original ixyz incident direction (y-like) + IF(icoor == 2)THEN ! z-like direction + nu = 10 + nv = 11 + nw = 9 + iratio = 12 + ELSEIF(icoor == 3)THEN ! x-like direction + nu = 2 + nv = 3 + nw = 1 + iratio = 4 + END IF + ELSEIF(ixyz == 3)THEN ! z-like direction + mu = 4 + mv = 2 + mw = 3 + ! Choose the coordinate direction to solve the Riemann problem + ! Permute nu,nv,nw for the original ixyz incident direction (z-like) + IF(icoor == 2)THEN ! x-like direction + nu = 3 + nv = 1 + nw = 2 + iratio = 4 + ELSEIF(icoor == 3)THEN ! y-like direction + nu = 7 + nv = 5 + nw = 6 + iratio = 8 + END IF + ENDIF + + ! Note that notation for u,v, and w reflects assumption that the + ! Riemann problems are in the x-direction with u in the normal + ! direction and v and w in the orthogonal directions, but with the + ! above definitions of mu, mv, and mw the routine also works with + ! ixyz=2 and ixyz = 3 + ! and returns, for example, f0 as the Godunov flux g0 for the + ! Riemann problems u_t + g(u)_y = 0 in the y-direction. + ! + ! Compute the Roe-averaged variables needed in the Roe solver. + + ! Loop over grid cell interfaces + DO i = 2-mbc, mx+mbc + IF (qr(1,i-1) <= 0.d0 .OR. ql(1,i) <= 0.d0) THEN + WRITE(*,*) i, mu, mv, mw + WRITE(*,'(5e12.4)') (qr(j,i-1),j=1,5) + WRITE(*,'(5e12.4)') (ql(j,i),j=1,5) + IF (ixyz == 1) WRITE(6,*) '*** (rpt3) rho <= 0 in x-sweep at ',i + IF (ixyz == 2) WRITE(6,*) '*** (rpt3) rho <= 0 in y-sweep at ',i + IF (ixyz == 3) WRITE(6,*) '*** (rpt3) rho <= 0 in z-sweep at ',i + WRITE(6,*) 'stopped with rho <= 0...' + STOP + ENDIF + rhsqrtl = SQRT(qr(1,i-1)) + rhsqrtr = SQRT(ql(1,i)) + pl = gamma1*(qr(5,i-1) - 0.5d0*(qr(mu,i-1)**2 + & + qr(mv,i-1)**2 + qr(mw,i-1)**2)/qr(1,i-1)) + pr = gamma1*(ql(5,i) - 0.5d0*(ql(mu,i)**2 + & + ql(mv,i)**2 + ql(mw,i)**2)/ql(1,i)) + rhsq2 = rhsqrtl + rhsqrtr + u(i) = (qr(mu,i-1)/rhsqrtl + ql(mu,i)/rhsqrtr) / rhsq2 + v(i) = (qr(mv,i-1)/rhsqrtl + ql(mv,i)/rhsqrtr) / rhsq2 + w(i) = (qr(mw,i-1)/rhsqrtl + ql(mw,i)/rhsqrtr) / rhsq2 + enth(i) = (((qr(5,i-1)+pl)/rhsqrtl & + + (ql(5,i)+pr)/rhsqrtr)) / rhsq2 + u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2 + asqrd = gamma1*(enth(i) - .5d0*u2v2w2(i)) + + IF (i>=0 .AND. i<=mx .AND. asqrd <= 0.d0) THEN + IF (ixyz == 1) WRITE(6,*) '*** (rpt3) a**2 <= 0 in x-sweep at ',i + IF (ixyz == 2) WRITE(6,*) '*** (rpt3) a**2 <= 0 in y-sweep at ',i + IF (ixyz == 3) WRITE(6,*) '*** (rpt3) a**2 <= 0 in z-sweep at ',i + WRITE(6,*) 'stopped with a**2 < 0...' + STOP + ENDIF + a(i) = SQRT(asqrd) + g1a2(i) = gamma1 / asqrd + euv(i) = enth(i) - u2v2w2(i) + END DO + + ! Solve the Riemann Problem + DO i = 2-mbc, mx+mbc + + nx = aux2(nu,i,2) + ny = aux2(nv,i,2) + nz = aux2(nw,i,2) + vn = u(i)*nx + v(i)*ny + w(i)*nz + IF(icoor == 2)THEN + areaRatioMinus = aux2(iratio,i-2+imp,2) + areaRatioPlus = aux3(iratio,i-2+imp,2) + ELSEIF(icoor == 3)THEN + areaRatioMinus = aux2(iratio,i-2+imp,2) + areaRatioPlus = aux2(iratio,i-2+imp,3) + ENDIF + ! Pick the largest normal to avoid using a singular direction + IF(ABS(nx) >= ABS(ny) .AND. ABS(nx) >= ABS(nz))THEN + nonsingularNormal = "x" + ELSEIF(ABS(ny) >= ABS(nx) .AND. ABS(ny) >= ABS(nz))THEN + nonsingularNormal = "y" + ELSEIF(ABS(nz) >= ABS(nx) .AND. ABS(nz) >= ABS(ny))THEN + nonsingularNormal = "z" + ELSE + WRITE(*,*)"Invalid Normal Vector..." + WRITE(*,*)"nx: ",nx,"ny: ",ny,"nz: ",nz + STOP + END IF + + beta(4) = g1a2(i)*(euv(i)*asdq(1,i) & + + u(i)*asdq(mu,i) + v(i)*asdq(mv,i) + w(i)*asdq(mw,i) - asdq(5,i)) + beta(5) = ( (a(i)-vn)*asdq(1,i) - a(i)*beta(4) & + + nx*asdq(mu,i) + ny*asdq(mv,i) + nz*asdq(mw,i) )/(2.d0*a(i)) + beta(1) = asdq(1,i) - beta(4) - beta(5) + + SELECT CASE(TRIM(nonsingularNormal)) + CASE("x") + beta(2) = ( (v(i)-vn*ny)*asdq(1,i) + nx*ny*asdq(mu,i) & + + (ny**2-1)*asdq(mv,i) + ny*nz*asdq(mw,i) )/nx + beta(3) = (-(w(i)-vn*nz)*asdq(1,i) - nx*nz*asdq(mu,i) & + - ny*nz*asdq(mv,i) - (nz**2-1)*asdq(mw,i) )/nx + CASE("y") + beta(2) = (-(u(i)-vn*nx)*asdq(1,i) - (nx**2-1)*asdq(mu,i) & + - ny*nx*asdq(mv,i) - nx*nz*asdq(mw,i) )/ny + beta(3) = ( (w(i)-vn*nz)*asdq(1,i) + nx*nz*asdq(mu,i) & + + ny*nz*asdq(mv,i) + (nz**2-1)*asdq(mw,i) )/ny + CASE("z") + beta(2) = ( (u(i)-vn*nx)*asdq(1,i) + (nx**2-1)*asdq(mu,i) & + + nx*ny*asdq(mv,i) + nz*nx*asdq(mw,i) )/nz + beta(3) = (-(v(i)-vn*ny)*asdq(1,i) - nx*ny*asdq(mu,i) & + - (ny**2-1)*asdq(mv,i) - nz*ny*asdq(mw,i) )/nz + END SELECT + + waveb(1,1) = beta(1) + waveb(mu,1) = beta(1)*(u(i)-a(i)*nx) + waveb(mv,1) = beta(1)*(v(i)-a(i)*ny) + waveb(mw,1) = beta(1)*(w(i)-a(i)*nz) + waveb(5,1) = beta(1)*(enth(i) - vn*a(i)) + sb(1) = vn - a(i) + + waveb(1,2) = beta(4) + SELECT CASE(TRIM(nonsingularNormal)) + CASE("x") + waveb(mu,2) = beta(4)*u(i) + beta(2)*ny - beta(3)*nz + waveb(mv,2) = beta(4)*v(i) - beta(2)*nx + waveb(mw,2) = beta(4)*w(i) + beta(3)*nx + waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & + + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(w(i)*nx-u(i)*nz) + CASE("y") + waveb(mu,2) = beta(4)*u(i) + beta(2)*ny + waveb(mv,2) = beta(4)*v(i) - beta(2)*nx + beta(3)*nz + waveb(mw,2) = beta(4)*w(i) - beta(3)*ny + waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & + + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(v(i)*nz-w(i)*ny) + CASE("z") + waveb(mu,2) = beta(4)*u(i) - beta(2)*nz + waveb(mv,2) = beta(4)*v(i) + beta(3)*nz + waveb(mw,2) = beta(4)*w(i) + beta(2)*nx - beta(3)*ny + waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & + + beta(2)*(w(i)*nx-u(i)*nz) + beta(3)*(v(i)*nz-w(i)*ny) + END SELECT + sb(2) = vn + + waveb(1,3) = beta(5) + waveb(mu,3) = beta(5)*(u(i) + a(i)*nx) + waveb(mv,3) = beta(5)*(v(i) + a(i)*ny) + waveb(mw,3) = beta(5)*(w(i) + a(i)*nz) + waveb(5,3) = beta(5)*(enth(i) + vn*a(i)) + sb(3) = vn + a(i) + + bmasdq(:,i) = 0.d0 + bpasdq(:,i) = 0.d0 + DO mws = 1,mwaves + IF (sb(mws) < 0.d0)THEN + waveb(:,mws) = waveb(:,mws)*areaRatioMinus + bmasdq(:,i) = bmasdq(:,i) + sb(mws)*waveb(:,mws) + ELSE + waveb(:,mws) = waveb(:,mws)*areaRatioPlus + bpasdq(:,i) = bpasdq(:,i) + sb(mws)*waveb(:,mws) + END IF + END DO + END DO + +END SUBROUTINE rpt3 diff --git a/examples/euler_gravity_3d/rptt3_euler_mapgrid.f90 b/examples/euler_gravity_3d/rptt3_euler_mapgrid.f90 new file mode 100644 index 000000000..9a4392769 --- /dev/null +++ b/examples/euler_gravity_3d/rptt3_euler_mapgrid.f90 @@ -0,0 +1,306 @@ +!----------------------------------------------------------------------------------- +! Name: rptt3(ixyz,icoor,imp,impt,maxm,meqn,mwaves,maux,mbc,mx,ql,qr, +! aux1,aux2,aux3,asdq,cmbsasdq,cpbsasdq) +! +! Description: Roe f-wave solver in the transverse direction for the +! 3D Euler equations in general geometries +! +! Inputs: +! ixyz : direction to take slice x-direction if ixyz=1 +! y-direction if ixyz=2. +! z-direction if ixyz=3. +! icoor : direction in which the transverse solve should be +! performed. icoor=2 (split in y-like direction) +! icoor=3 (split in z-like direction) +! ixyz=1, icoor=3 means bsasdq=B*A*\Dq, split in z into: +! cmbsasdq = C^-B*A*\Dq, +! cpbsasdq = C^+B*A*\Dq. +! ixyz=2, icoor=3 means bsasdq=C*B*\Dq, split in x into: +! cmbsasdq = A^-C*B*\Dq, +! cpbsasdq = A^+C*B*\Dq. +! imp : index for aux arrays +! impt : index for aux arrays +! maxm : max number of grid cells (less the ghost cells) +! meqn : number of equations in the system +! mwaves : nuber of waves in the system +! maux : number of auxilary equations +! mbc : number of ghost cells on either end +! mx : number of elements +! ql : state vector at left edge of each cell +! Note the i'th Riemann problem has left state qr(:,i-1) +! qr : state vector at right edge of each cell +! Note the i'th Riemann problem has right state ql(:,i) +! aux1 : +! aux2 : +! aux3 : +! asdq : array of flux differences (A*\Dq). asdq(:,i) is the +! flux difference propagating away from the interface +! between cells i-1 and i. Note: asdq represents +! B*\Dq if ixyz=2 or C*\Dq if ixyz=3. +! +! Outputs: +! cmbsasdq : left-going flux differences +! cpbsasdq : right-going flux differences +! +! Adapted from rptt3_euler.f90 in $CLAWHOME/riemann/src +!----------------------------------------------------------------------------------- +SUBROUTINE rptt3(ixyz,icoor,imp,impt,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,& + aux1,aux2,aux3,bsasdq,cmbsasdq,cpbsasdq) + + IMPLICIT NONE + + ! Input + INTEGER, INTENT(IN) :: ixyz,icoor,imp,impt,maxm,meqn,mwaves,maux,mbc,mx + REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(IN) :: ql,qr,bsasdq + REAL(kind=8), DIMENSION(maux,1-mbc:maxm+mbc,3), INTENT(IN) :: aux1,aux2,aux3 + + ! Output + REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(INOUT) :: cmbsasdq,cpbsasdq + + ! Local Storage + REAL(kind=8), DIMENSION(meqn) :: beta + REAL(kind=8), DIMENSION(-1:maxm+mbc) :: u2v2w2,u,v,w,enth,a,g1a2,euv + REAL(kind=8) :: gamma1,asqrd,pl,pr,rhsq2,rhsqrtl,rhsqrtr,nx,ny,nz,vn,& + areaRatioMinus,areaRatioPlus + REAL(kind=8), DIMENSION(meqn,mwaves) :: waveb + REAL(kind=8), DIMENSION(mwaves) :: sb + INTEGER :: i,j,mu=2,mv=3,mw=4,mws,nu=1,nv=2,nw=3,iratio=4 + CHARACTER(20) :: nonsingularNormal + + ! Common block storage for ideal gas constant + REAL(kind=8) :: gamma + COMMON /cparam/ gamma + + ! Set (gamma-1) + gamma1 = gamma - 1.d0 + + ! Make sure the number of waves is 3 for the 3D Euler equations + IF (mwaves /= 3) THEN + WRITE(6,*) '*** Should have mwaves=3 for this Riemann solver' + STOP + ENDIF + + ! Set mu to point to the component of the system that corresponds to momentum in + ! the direction of this slice, mv and mw to the orthogonal momentums. + ! + ! ixyz indicates the direction of the original Riemann solve, + ! called the x-like direction in the table below: + ! + ! x-like direction y-like direction z-like direction + ! ixyz=1: x y z + ! ixyz=2: y z x + ! ixyz=3: z x y + IF(ixyz == 1)THEN ! x-like direction + mu = 2 + mv = 3 + mw = 4 + ! Choose the coordinate direction to solve the Riemann problem + ! Permute nu,nv,nw for the original ixyz incident direction (x-like) + IF(icoor == 2)THEN ! y-like direction + nu = 5 + nv = 6 + nw = 7 + iratio = 8 + ELSEIF(icoor == 3)THEN ! z-like direction + nu = 9 + nv = 10 + nw = 11 + iratio = 12 + END IF + ELSEIF(ixyz == 2)THEN ! y-like direction + mu = 3 + mv = 4 + mw = 2 + ! Choose the coordinate direction to solve the Riemann problem + ! Permute nu,nv,nw for the original ixyz incident direction (y-like) + IF(icoor == 2)THEN ! z-like direction + nu = 10 + nv = 11 + nw = 9 + iratio = 12 + ELSEIF(icoor == 3)THEN ! x-like direction + nu = 2 + nv = 3 + nw = 1 + iratio = 4 + END IF + ELSEIF(ixyz == 3)THEN ! z-like direction + mu = 4 + mv = 2 + mw = 3 + ! Choose the coordinate direction to solve the Riemann problem + ! Permute nu,nv,nw for the original ixyz incident direction (z-like) + IF(icoor == 2)THEN ! x-like direction + nu = 3 + nv = 1 + nw = 2 + iratio = 4 + ELSEIF(icoor == 3)THEN ! y-like direction + nu = 7 + nv = 5 + nw = 6 + iratio = 8 + END IF + ENDIF + + ! Note that notation for u,v, and w reflects assumption that the + ! Riemann problems are in the x-direction with u in the normal + ! direction and v and w in the orthogonal directions, but with the + ! above definitions of mu, mv, and mw the routine also works with + ! ixyz=2 and ixyz = 3 + ! and returns, for example, f0 as the Godunov flux g0 for the + ! Riemann problems u_t + g(u)_y = 0 in the y-direction. + ! + ! Compute the Roe-averaged variables needed in the Roe solver. + + ! Loop over grid cell interfaces + DO i = 2-mbc, mx+mbc + IF (qr(1,i-1) <= 0.d0 .OR. ql(1,i) <= 0.d0) THEN + WRITE(*,*) i, mu, mv, mw + WRITE(*,'(5e12.4)') (qr(j,i-1),j=1,5) + WRITE(*,'(5e12.4)') (ql(j,i),j=1,5) + IF (ixyz == 1) WRITE(6,*) '*** (rptt3) rho <= 0 in x-sweep at ',i + IF (ixyz == 2) WRITE(6,*) '*** (rptt3) rho <= 0 in y-sweep at ',i + IF (ixyz == 3) WRITE(6,*) '*** (rptt3) rho <= 0 in z-sweep at ',i + WRITE(6,*) 'stopped with rho <= 0...' + STOP + ENDIF + rhsqrtl = SQRT(qr(1,i-1)) + rhsqrtr = SQRT(ql(1,i)) + pl = gamma1*(qr(5,i-1) - 0.5d0*(qr(mu,i-1)**2 + & + qr(mv,i-1)**2 + qr(mw,i-1)**2)/qr(1,i-1)) + pr = gamma1*(ql(5,i) - 0.5d0*(ql(mu,i)**2 + & + ql(mv,i)**2 + ql(mw,i)**2)/ql(1,i)) + rhsq2 = rhsqrtl + rhsqrtr + u(i) = (qr(mu,i-1)/rhsqrtl + ql(mu,i)/rhsqrtr) / rhsq2 + v(i) = (qr(mv,i-1)/rhsqrtl + ql(mv,i)/rhsqrtr) / rhsq2 + w(i) = (qr(mw,i-1)/rhsqrtl + ql(mw,i)/rhsqrtr) / rhsq2 + enth(i) = (((qr(5,i-1)+pl)/rhsqrtl & + + (ql(5,i)+pr)/rhsqrtr)) / rhsq2 + u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2 + asqrd = gamma1*(enth(i) - .5d0*u2v2w2(i)) + + IF (i>=0 .AND. i<=mx .AND. asqrd <= 0.d0) THEN + IF (ixyz == 1) WRITE(6,*) '*** (rptt3) a**2 <= 0 in x-sweep at ',i + IF (ixyz == 2) WRITE(6,*) '*** (rptt3) a**2 <= 0 in y-sweep at ',i + IF (ixyz == 3) WRITE(6,*) '*** (rptt3) a**2 <= 0 in z-sweep at ',i + WRITE(6,*) 'stopped with a**2 < 0...' + STOP + ENDIF + a(i) = SQRT(asqrd) + g1a2(i) = gamma1 / asqrd + euv(i) = enth(i) - u2v2w2(i) + END DO + + ! Solve the Riemann Problem + DO i = 2-mbc, mx+mbc + + nx = aux2(nu,i,2) + ny = aux2(nv,i,2) + nz = aux2(nw,i,2) + vn = u(i)*nx + v(i)*ny + w(i)*nz + IF(icoor == 2)THEN + IF(impt==1)THEN + areaRatioMinus = aux2(iratio,i-2+imp,1) + areaRatioPlus = aux3(iratio,i-2+imp,1) + ELSEIF(impt==2)THEN + areaRatioMinus = aux2(iratio,i-2+imp,3) + areaRatioPlus = aux3(iratio,i-2+imp,3) + END IF + ELSEIF(icoor == 3)THEN + IF(impt == 1)THEN + areaRatioMinus = aux1(iratio,i-2+imp,2) + areaRatioPlus = aux1(iratio,i-2+imp,3) + ELSEIF(impt == 2)THEN + areaRatioMinus = aux3(iratio,i-2+imp,2) + areaRatioPlus = aux3(iratio,i-2+imp,3) + END IF + END IF + ! Pick the largest normal to avoid using a singular direction + IF(ABS(nx) >= ABS(ny) .AND. ABS(nx) >= ABS(nz))THEN + nonsingularNormal = "x" + ELSEIF(ABS(ny) >= ABS(nx) .AND. ABS(ny) >= ABS(nz))THEN + nonsingularNormal = "y" + ELSEIF(ABS(nz) >= ABS(nx) .AND. ABS(nz) >= ABS(ny))THEN + nonsingularNormal = "z" + ELSE + WRITE(*,*)"Invalid Normal Vector..." + WRITE(*,*)"nx: ",nx,"ny: ",ny,"nz: ",nz + STOP + END IF + + beta(4) = g1a2(i)*(euv(i)*bsasdq(1,i) & + + u(i)*bsasdq(mu,i) + v(i)*bsasdq(mv,i) + w(i)*bsasdq(mw,i) - bsasdq(5,i)) + beta(5) = ( (a(i)-vn)*bsasdq(1,i) - a(i)*beta(4) & + + nx*bsasdq(mu,i) + ny*bsasdq(mv,i) + nz*bsasdq(mw,i) )/(2.d0*a(i)) + beta(1) = bsasdq(1,i) - beta(4) - beta(5) + + SELECT CASE(TRIM(nonsingularNormal)) + CASE("x") + beta(2) = ( (v(i)-vn*ny)*bsasdq(1,i) + nx*ny*bsasdq(mu,i) & + + (ny**2-1)*bsasdq(mv,i) + ny*nz*bsasdq(mw,i) )/nx + beta(3) = (-(w(i)-vn*nz)*bsasdq(1,i) - nx*nz*bsasdq(mu,i) & + - ny*nz*bsasdq(mv,i) - (nz**2-1)*bsasdq(mw,i) )/nx + CASE("y") + beta(2) = (-(u(i)-vn*nx)*bsasdq(1,i) - (nx**2-1)*bsasdq(mu,i) & + - ny*nx*bsasdq(mv,i) - nx*nz*bsasdq(mw,i) )/ny + beta(3) = ( (w(i)-vn*nz)*bsasdq(1,i) + nx*nz*bsasdq(mu,i) & + + ny*nz*bsasdq(mv,i) + (nz**2-1)*bsasdq(mw,i) )/ny + CASE("z") + beta(2) = ( (u(i)-vn*nx)*bsasdq(1,i) + (nx**2-1)*bsasdq(mu,i) & + + nx*ny*bsasdq(mv,i) + nz*nx*bsasdq(mw,i) )/nz + beta(3) = (-(v(i)-vn*ny)*bsasdq(1,i) - nx*ny*bsasdq(mu,i) & + - (ny**2-1)*bsasdq(mv,i) - nz*ny*bsasdq(mw,i) )/nz + END SELECT + + waveb(1,1) = beta(1) + waveb(mu,1) = beta(1)*(u(i)-a(i)*nx) + waveb(mv,1) = beta(1)*(v(i)-a(i)*ny) + waveb(mw,1) = beta(1)*(w(i)-a(i)*nz) + waveb(5,1) = beta(1)*(enth(i) - vn*a(i)) + sb(1) = vn - a(i) + + waveb(1,2) = beta(4) + SELECT CASE(TRIM(nonsingularNormal)) + CASE("x") + waveb(mu,2) = beta(4)*u(i) + beta(2)*ny - beta(3)*nz + waveb(mv,2) = beta(4)*v(i) - beta(2)*nx + waveb(mw,2) = beta(4)*w(i) + beta(3)*nx + waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & + + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(w(i)*nx-u(i)*nz) + CASE("y") + waveb(mu,2) = beta(4)*u(i) + beta(2)*ny + waveb(mv,2) = beta(4)*v(i) - beta(2)*nx + beta(3)*nz + waveb(mw,2) = beta(4)*w(i) - beta(3)*ny + waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & + + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(v(i)*nz-w(i)*ny) + CASE("z") + waveb(mu,2) = beta(4)*u(i) - beta(2)*nz + waveb(mv,2) = beta(4)*v(i) + beta(3)*nz + waveb(mw,2) = beta(4)*w(i) + beta(2)*nx - beta(3)*ny + waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & + + beta(2)*(w(i)*nx-u(i)*nz) + beta(3)*(v(i)*nz-w(i)*ny) + END SELECT + sb(2) = vn + + waveb(1,3) = beta(5) + waveb(mu,3) = beta(5)*(u(i) + a(i)*nx) + waveb(mv,3) = beta(5)*(v(i) + a(i)*ny) + waveb(mw,3) = beta(5)*(w(i) + a(i)*nz) + waveb(5,3) = beta(5)*(enth(i) + vn*a(i)) + sb(3) = vn + a(i) + + cmbsasdq(:,i) = 0.d0 + cpbsasdq(:,i) = 0.d0 + DO mws = 1,mwaves + IF (sb(mws) < 0.d0)THEN + waveb(:,mws) = waveb(:,mws)*areaRatioMinus + cmbsasdq(:,i) = cmbsasdq(:,i) + sb(mws)*waveb(:,mws) + ELSE + waveb(:,mws) = waveb(:,mws)*areaRatioPlus + cpbsasdq(:,i) = cpbsasdq(:,i) + sb(mws)*waveb(:,mws) + END IF + END DO + END DO + +END SUBROUTINE rptt3 From 29edc6a7ed277b95c1088a10647a87776d2ed55f Mon Sep 17 00:00:00 2001 From: Carlos Munoz Moncayo Date: Thu, 4 Jul 2024 12:33:48 +0300 Subject: [PATCH 2/7] Update meson.build --- src/pyclaw/meson.build | 47 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/src/pyclaw/meson.build b/src/pyclaw/meson.build index f99edd2e8..3f7f3c038 100644 --- a/src/pyclaw/meson.build +++ b/src/pyclaw/meson.build @@ -147,6 +147,7 @@ examples = { 'test_sedov_and_hdf.py', ], 'euler_gravity_3d' : [ + '__init__.py', 'plotCreateVisitXDMF.py', 'rising_hot_sphere.py', 'rising_hot_sphere_spherical.py', @@ -325,3 +326,49 @@ py.extension_module( subdir: 'clawpack/pyclaw/examples/advection_reaction_2d', install : true ) + +# 3D Euler with gravity + +ext_name = 'euler_3d_gmap' +ext_srcs = [ + '../../examples/euler_gravity_3d/rpn3_euler_mapgrid.f90', + '../../examples/euler_gravity_3d/rpt3_euler_mapgrid.f90', + '../../examples/euler_gravity_3d/rptt3_euler_mapgrid.f90', +] +f2py_srcs = custom_target( + 'f2py_3d_gmap', + command: [f2py, ext_name], + input: ext_srcs, + output: [ext_name + 'module.c', ext_name + '-f2pywrappers.f'], + build_by_default: true, +) + +py.extension_module( + ext_name, [ext_srcs, f2py_srcs], + incdir_f2py / 'fortranobject.c', + include_directories: inc_np, + dependencies : py_dep, + subdir: 'clawpack/pyclaw/examples/euler_gravity_3d', + install : true +) + +ext_name = 'mappedGrid' +ext_srcs = [ + '../../examples/euler_gravity_3d/euler3d_mappedGrid.f90', +] +f2py_srcs = custom_target( + 'f2py_mappedGrid', + command: [f2py, ext_name], + input: ext_srcs, + output: [ext_name + 'module.c', ext_name + '-f2pywrappers2.f90'], + build_by_default: true, +) + +py.extension_module( + ext_name, [ext_srcs, f2py_srcs], + incdir_f2py / 'fortranobject.c', + include_directories: inc_np, + dependencies : py_dep, + subdir: 'clawpack/pyclaw/examples/euler_gravity_3d', + install : true +) \ No newline at end of file From 69090d0ddd6dec337e96e8955e2840541ebe8856 Mon Sep 17 00:00:00 2001 From: Carlos Munoz Moncayo Date: Thu, 4 Jul 2024 12:36:05 +0300 Subject: [PATCH 3/7] Clean old comments --- examples/euler_gravity_3d/rising_hot_sphere.py | 2 -- examples/euler_gravity_3d/rising_hot_sphere_spherical.py | 2 -- 2 files changed, 4 deletions(-) diff --git a/examples/euler_gravity_3d/rising_hot_sphere.py b/examples/euler_gravity_3d/rising_hot_sphere.py index 24e54a330..b95d006f6 100755 --- a/examples/euler_gravity_3d/rising_hot_sphere.py +++ b/examples/euler_gravity_3d/rising_hot_sphere.py @@ -12,7 +12,6 @@ """ import numpy as np from clawpack.pyclaw.examples.euler_gravity_3d.mappedGrid import euler3d_mappedgrid as mg -# from .mappedGrid import euler3d_mappedgrid as mg try: from mpi4py import MPI @@ -178,7 +177,6 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ import logging solver.logger.setLevel(logging.DEBUG) from clawpack.pyclaw.examples.euler_gravity_3d import euler_3d_gmap - # from . import euler_3d_gmap solver.rp = euler_3d_gmap solver.num_eqn = 5 solver.num_waves = 3 diff --git a/examples/euler_gravity_3d/rising_hot_sphere_spherical.py b/examples/euler_gravity_3d/rising_hot_sphere_spherical.py index 41ce20c08..58e62a335 100755 --- a/examples/euler_gravity_3d/rising_hot_sphere_spherical.py +++ b/examples/euler_gravity_3d/rising_hot_sphere_spherical.py @@ -12,7 +12,6 @@ """ import numpy as np from clawpack.pyclaw.examples.euler_gravity_3d.mappedGrid import euler3d_mappedgrid as mg -# from .mappedGrid import euler3d_mappedgrid as mg # Test for MPI, and set sizes accordingly try: @@ -175,7 +174,6 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ import logging solver.logger.setLevel(logging.DEBUG) from clawpack.pyclaw.examples.euler_gravity_3d import euler_3d_gmap - #from . import euler_3d_gmap solver.rp = euler_3d_gmap solver.num_eqn = 5 solver.num_waves = 3 From 2b815b93bf14db64e013e533e52f16cf05c9ad00 Mon Sep 17 00:00:00 2001 From: Carlos Munoz Moncayo Date: Thu, 4 Jul 2024 16:43:20 +0300 Subject: [PATCH 4/7] Leave source in the Riemann repository --- examples/euler_gravity_3d/Makefile | 5 +- .../euler_gravity_3d/euler3d_mappedGrid.f90 | 679 ------------------ .../euler_gravity_3d/rpn3_euler_mapgrid.f90 | 312 -------- .../euler_gravity_3d/rpt3_euler_mapgrid.f90 | 295 -------- .../euler_gravity_3d/rptt3_euler_mapgrid.f90 | 306 -------- 5 files changed, 3 insertions(+), 1594 deletions(-) delete mode 100644 examples/euler_gravity_3d/euler3d_mappedGrid.f90 delete mode 100644 examples/euler_gravity_3d/rpn3_euler_mapgrid.f90 delete mode 100644 examples/euler_gravity_3d/rpt3_euler_mapgrid.f90 delete mode 100644 examples/euler_gravity_3d/rptt3_euler_mapgrid.f90 diff --git a/examples/euler_gravity_3d/Makefile b/examples/euler_gravity_3d/Makefile index 1ca2c7a93..182410491 100644 --- a/examples/euler_gravity_3d/Makefile +++ b/examples/euler_gravity_3d/Makefile @@ -1,14 +1,15 @@ F2PY = f2py +RIEMANN_SOURCE = ../../../riemann/src/ all: $(MAKE) euler_3d_gmap.so $(MAKE) mappedGrid.so -euler_3d_gmap.so: rpn3_euler_mapgrid.f90 rpt3_euler_mapgrid.f90 rptt3_euler_mapgrid.f90 +euler_3d_gmap.so: $(RIEMANN_SOURCE)/rpn3_euler_mapgrid.f90 $(RIEMANN_SOURCE)/rpt3_euler_mapgrid.f90 $(RIEMANN_SOURCE)/rptt3_euler_mapgrid.f90 f2py -m $(basename $(notdir $@)) -c $^ -mappedGrid.so: euler3d_mappedGrid.f90 +mappedGrid.so: $(RIEMANN_SOURCE)/euler3d_mappedGrid.f90 f2py -m $(basename $(notdir $@)) -c $^ clean: diff --git a/examples/euler_gravity_3d/euler3d_mappedGrid.f90 b/examples/euler_gravity_3d/euler3d_mappedGrid.f90 deleted file mode 100644 index 55d149fdf..000000000 --- a/examples/euler_gravity_3d/euler3d_mappedGrid.f90 +++ /dev/null @@ -1,679 +0,0 @@ -MODULE euler3d_mappedgrid - - REAL(kind=8), PARAMETER :: pi=3.1415926535897932385,twopi=2*pi - REAL(kind=8), PARAMETER :: rEarth=637120000.0 - -CONTAINS - !------------------------------------------------------------------------- - ! Name: - ! volumeOfHex(x,y,z): - ! - ! Description: - ! Compute Volume of an LD Hexahedron given the 8 vertices - ! from paper: - ! @book{Grandy_1997, - ! title={Efficient computation of volume of hexahedral cells}, - ! url={http://www.osti.gov/scitech/servlets/purl/632793}, - ! DOI={10.2172/632793}, - ! abstractNote={This report describes an efficient method to compute - ! the volume of hexahedral cells used in three-dimensional - ! hydrodynamics simulation. Two common methods for - ! creating the hexahedron using triangular boundaries are - ! considered.}, - ! author={Grandy, J.}, year={1997}, month={Oct}} - ! https://doi.org/10.2172/632793 - ! - ! 6*V_LD = - ! |x6-x0 x3-x0 x2-x7| |x6-x0 x4-x0 x7-x5| |x6-x0 x1-x0 x5-x2| - ! |y6-y0 y3-y0 y2-y7| + |y6-y0 y4-y0 y7-y5| + |y6-y0 y1-y0 y5-y2| - ! |z6-z0 z3-z0 z2-z7| |z6-z0 z4-z0 z7-z5| |z6-z0 z1-z0 z5-z2| - ! - ! Inputs: - ! x[:],y[:],z[:] : x,y, and z coordinates of eight hexahedron vertices - ! - ! Output: - ! volumeOfHex : Volume of LD hexahedron - !------------------------------------------------------------------------- - REAL(kind=8) FUNCTION volumeOfHex(x,y,z) - - IMPLICIT NONE - - REAL(kind=8), DIMENSION(8), INTENT(IN) :: x,y,z - - ! Compute Volume of Hexahedron - volumeOfHex = & - (x(7)-x(1))*((y(4)-y(1))*(z(3)-z(8)) - (y(3)-y(8))*(z(4)-z(1)))+& - (x(4)-x(1))*((y(3)-y(8))*(z(7)-z(1)) - (y(7)-y(1))*(z(3)-z(8)))+& - (x(3)-x(8))*((y(7)-y(1))*(z(4)-z(1)) - (y(4)-y(1))*(z(7)-z(1)))& - +& - (x(7)-x(1))*((y(5)-y(1))*(z(8)-z(6)) - (y(8)-y(6))*(z(5)-z(1)))+& - (x(5)-x(1))*((y(8)-y(6))*(z(7)-z(1)) - (y(7)-y(1))*(z(8)-z(6)))+& - (x(8)-x(6))*((y(7)-y(1))*(z(5)-z(1)) - (y(5)-y(1))*(z(7)-z(1)))& - +& - (x(7)-x(1))*((y(2)-y(1))*(z(6)-z(3)) - (y(6)-y(3))*(z(2)-z(1)))+& - (x(2)-x(1))*((y(6)-y(3))*(z(7)-z(1)) - (y(7)-y(1))*(z(6)-z(3)))+& - (x(6)-x(3))*((y(7)-y(1))*(z(2)-z(1)) - (y(2)-y(1))*(z(7)-z(1))) - volumeOfHex = ABS(volumeOfHex)/6.d0 - RETURN - END FUNCTION volumeOfHex - - !--------------------------------------------------------------------------- - ! Name: mapc2pWrapper(xc,yc,zc,mx,my,mz,xyz0,xyzN,mapType,xp,yp,zp) - ! - ! Description: mapping from computational to physical space - ! - ! Inputs: - ! xc,yc,zc : computational coordinate values - ! mx,my,mz : grid sizes - ! xyz0,xyzN : altitude at bottom and top of grid - ! mapType : type of mapping to compute - ! - ! Outputs: - ! xp,yp,zp : physical coordinate values - !--------------------------------------------------------------------------- - SUBROUTINE mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc - REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN - INTEGER, INTENT(IN) :: mz - CHARACTER(20) :: mapType - - ! Output - REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp - REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp - REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp - - ! Compute Physical Grid From Computational - SELECT CASE(mapType) - CASE("Identity") - CALL mapc2pIdentity(xc,yc,zc,xp,yp,zp) - CASE("ZeroToOne") - CALL mapc2pZeroToOne(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - CASE("Rotate") - CALL mapc2pRotate(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - CASE("Spherical") - CALL mapc2pSpherical(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - CASE DEFAULT - WRITE(*,*)"Invalid mapType: ",mapType - STOP - END SELECT - - RETURN - END SUBROUTINE mapc2pWrapper - - !--------------------------------------------------------------------------- - ! Name: mapc2pIdentity(xc,yc,zc,xp,yp,zp) - ! - ! Description: Identity mapping from computational to physical space - ! - ! Inputs: - ! xc,yc,zc : computational coordinate values - ! Outputs: - ! xp,yp,zp : physical coordinate values - !--------------------------------------------------------------------------- - SUBROUTINE mapc2pIdentity(xc,yc,zc,xp,yp,zp) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc - - ! Output - REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp - REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp - REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp - - ! x,xi coordinate - xp = xc - - ! y,eta coordinate - yp = yc - - ! z,phi coordinate - zp = zc - - RETURN - END SUBROUTINE mapc2pIdentity - - !--------------------------------------------------------------------------- - ! Name: mapc2pZeroToOne(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - ! - ! Description: mapping from computational to physical space using a - ! computational grid from 0 to 1 - ! - ! Inputs: - ! xc,yc,zc : computational coordinate values - ! xyz0,xyzN : physical coordinates at bottom and top of grid - ! - ! Outputs: - ! xp,yp,zp : physical coordinate values - !--------------------------------------------------------------------------- - SUBROUTINE mapc2pZeroToOne(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc - REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN - - ! Output - REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp - REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp - REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp - - ! x,xi coordinate - xp = xyz0(1) + xc*(xyzN(1)-xyz0(1)) - - ! y,eta coordinate - yp = xyz0(2) + yc*(xyzN(2)-xyz0(2)) - - ! z,phi coordinate - zp = xyz0(3) + zc*(xyzN(3)-xyz0(3)) - - RETURN - END SUBROUTINE mapc2pZeroToOne - - !--------------------------------------------------------------------------- - ! Name: mapc2pRotate(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - ! - ! Description: mapping from computational to physical space using a - ! computational grid from 0 to 1 and rotated by angle theta - ! - ! Inputs: - ! xc,yc,zc : computational coordinate values - ! xyz0,xyzN : physical coordinates at bottom and top of grid - ! - ! Outputs: - ! xp,yp,zp : physical coordinate values - !--------------------------------------------------------------------------- - SUBROUTINE mapc2pRotate(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc - REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN - - ! Output - REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp - REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp - REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp - - ! Local - REAL(kind=8) :: theta - REAL(kind=8), DIMENSION(SIZE(xc)) :: xp0 - REAL(kind=8), DIMENSION(SIZE(yc)) :: yp0 - REAL(kind=8), DIMENSION(SIZE(zc)) :: zp0 - - theta=0.d0*pi/180.d0 - - ! x,xi coordinate - xp = xyz0(1) + xc*(xyzN(1)-xyz0(1)) - - ! y,eta coordinate - yp = xyz0(2) + yc*(xyzN(2)-xyz0(2)) - - ! z,phi coordinate - zp = xyz0(3) + zc*(xyzN(3)-xyz0(3)) - - xp0 = xp - yp0 = yp - zp0 = zp - ! Rotate about x-axis by angle theta - xp = xp0 - yp = yp0*COS(theta) - zp0*SIN(theta) - zp = yp0*SIN(theta) + zp0*COS(theta) - - xp0 = xp - yp0 = yp - zp0 = zp - ! Rotate about y-axis by angle theta - xp = xp0*COS(theta) + zp0*SIN(theta) - yp = yp0 - zp =-xp0*SIN(theta) + zp0*COS(theta) - - xp0 = xp - yp0 = yp - zp0 = zp - ! Rotate about z-axis by angle theta - xp = xp0*COS(theta) - yp0*SIN(theta) - yp = xp0*SIN(theta) + yp0*COS(theta) - zp = zp0 - - RETURN - END SUBROUTINE mapc2pRotate - - !--------------------------------------------------------------------------- - ! Name: mapc2pSpherical(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - ! - ! Description: mapping from computational to physical space using a - ! computational grid from 0 to 1 and then using a - ! spherical transformation - ! - ! Inputs: - ! xc,yc,zc : computational coordinate values - ! xyz0,xyzN : physical coordinates at bottom and top of grid - ! - ! Outputs: - ! xp,yp,zp : physical coordinate values - !--------------------------------------------------------------------------- - SUBROUTINE mapc2pSpherical(xc,yc,zc,xyz0,xyzN,xp,yp,zp) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(:), INTENT(IN) :: xc,yc,zc - REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN - - ! Output - REAL(kind=8), DIMENSION(SIZE(xc)), INTENT(OUT) :: xp - REAL(kind=8), DIMENSION(SIZE(yc)), INTENT(OUT) :: yp - REAL(kind=8), DIMENSION(SIZE(zc)), INTENT(OUT) :: zp - - ! Local - REAL(kind=8), DIMENSION(SIZE(xc)) :: theta - REAL(kind=8), DIMENSION(SIZE(yc)) :: phi - REAL(kind=8), DIMENSION(SIZE(zc)) :: r - - ! x,xi coordinate - xp = xyz0(1) + xc*(xyzN(1)-xyz0(1)) - - ! y,eta coordinate - yp = xyz0(2) + yc*(xyzN(2)-xyz0(2)) - - ! z,phi coordinate - zp = xyz0(3) + zc*(xyzN(3)-xyz0(3)) - - ! Spherical to Cartesian Coordinates - theta = xp - phi = yp - r = zp - xp = r*SIN(theta)*COS(phi) - yp = r*SIN(theta)*SIN(phi) - zp = r*COS(theta) - - RETURN - END SUBROUTINE mapc2pSpherical - - !--------------------------------------------------------------------------- - ! Name: determinant(A,detA) - ! - ! Description: Compute the determinant of a 3x3 array - ! - ! Inputs: A : 3x3 Array - ! - ! Outputs: detA : determinant of A - !--------------------------------------------------------------------------- - SUBROUTINE determinant(A,detA) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(3,3), INTENT(IN) :: A - - ! Output - REAL(kind=8), INTENT(OUT) :: detA - - detA = A(1,1)*A(2,2)*A(3,3) & - + A(2,1)*A(3,2)*A(1,3) & - + A(3,1)*A(1,2)*A(2,3) & - - A(3,1)*A(2,2)*A(1,3) & - - A(2,1)*A(1,2)*A(3,3) & - - A(1,1)*A(3,2)*A(2,3) - - RETURN - END SUBROUTINE determinant - - !--------------------------------------------------------------------------- - ! Name: unitNormal(a,b,c,n) - ! - ! Description: Compute the unit normal of given 3 indices in 3D space - ! - ! Inputs: a,b,c : indices with x,y, and z coordinates - ! - ! Outputs: n : unit normals of plane defined from points a,b,c - !--------------------------------------------------------------------------- - SUBROUTINE unitNormal(a,b,c,n) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(3), INTENT(IN) :: a,b,c - - ! Output - REAL(kind=8), DIMENSION(3), INTENT(OUT) :: n - - ! Local - REAL(kind=8), DIMENSION(3,3) :: D - REAL(kind=8) :: magnitude - - D(:,1) = (/1.d0,a(2),a(3)/) - D(:,2) = (/1.d0,b(2),b(3)/) - D(:,3) = (/1.d0,c(2),c(3)/) - CALL determinant(D,n(1)) - D(:,1) = (/a(1),1.d0,a(3)/) - D(:,2) = (/b(1),1.d0,b(3)/) - D(:,3) = (/c(1),1.d0,c(3)/) - CALL determinant(D,n(2)) - D(:,1) = (/a(1),a(2),1.d0/) - D(:,2) = (/b(1),b(2),1.d0/) - D(:,3) = (/c(1),c(2),1.d0/) - CALL determinant(D,n(3)) - - magnitude = SQRT(n(1)**2 + n(2)**2 + n(3)**2) - n(1) = n(1)/magnitude - n(2) = n(2)/magnitude - n(3) = n(3)/magnitude - - RETURN - END SUBROUTINE unitNormal - - !--------------------------------------------------------------------------- - ! Name: dot(a,b,adotb) - ! - ! Description: Compute the dot product of vectors a and b - ! - ! Inputs: a,b : 3D vectors - ! - ! Outputs: adotb : dot product of a and b - !--------------------------------------------------------------------------- - SUBROUTINE dot(a,b,adotb) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(3), INTENT(IN) :: a,b - - ! Output - REAL(kind=8), INTENT(OUT) :: adotb - - adotb = a(1)*b(1) + a(2)*b(2) + a(3)*b(3) - - RETURN - END SUBROUTINE dot - - !--------------------------------------------------------------------------- - ! Name: cross(a,b,c) - ! - ! Description: Compute the cross product of vectors a and b - ! - ! Inputs: a,b : 3D vectors - ! - ! Outputs: c : cross product of a and b - !--------------------------------------------------------------------------- - SUBROUTINE cross(a,b,c) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(3), INTENT(IN) :: a,b - - ! Output - REAL(kind=8), DIMENSION(3), INTENT(OUT) :: c - - c(1) = a(2)*b(3) - a(3)*b(2) - c(2) = a(3)*b(1) - a(1)*b(3) - c(3) = a(1)*b(2) - a(2)*b(1) - - RETURN - END SUBROUTINE cross - - !--------------------------------------------------------------------------- - ! Name: areaPolygon(poly) - ! - ! Description: Compute the area of a polygon given its vertices - ! - ! Inputs: poly : vertices of polygon - ! - ! Outputs: area : area of polygon - !--------------------------------------------------------------------------- - SUBROUTINE areaPolygon(poly,polyArea) - - IMPLICIT NONE - - ! Input - REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: poly - - ! Output - REAL(kind=8), INTENT(OUT) :: polyArea - - ! Local - INTEGER :: i,np - REAL(kind=8), DIMENSION(3) :: total,prod,vi1,vi2,n - - ! Compute Area of a Polygon Given its Vertices - np = SIZE(poly,2) - IF(np < 3)THEN - WRITE(*,*)"This is not a polygon, no area" - STOP - END IF - total=0 - DO i = 1,np - vi1 = poly(:,i) - IF(i == np)THEN - vi2 = poly(:,1) - ELSE - vi2 = poly(:,i+1) - END IF - CALL cross(vi1,vi2,prod) - total(1) = total(1) + prod(1) - total(2) = total(2) + prod(2) - total(3) = total(3) + prod(3) - END DO - CALL unitNormal(poly(:,1),poly(:,2),poly(:,3),n) - CALL dot(total,n,polyArea) - polyArea = ABS(polyArea*0.5d0) - - RETURN - END SUBROUTINE areaPolygon - - !--------------------------------------------------------------------------- - ! Name: setAuxiliaryVariables(num_aux,mbc,mx,my,mz,xlower,ylower,zlower,& - ! dxc,dyc,dzc,xyz0,xyzN,aux) - ! - ! Description: Compute auxiliary variables for non-uniform/mapped grids - ! - ! Inputs: - ! xc,yc,zc : computational coordinate values - ! mx,my,mz : grid sizes - ! xyz0,xyzN : altitude at bottom and top of grid - ! - ! Outputs: - ! xp,yp,zp : physical coordinate values - !--------------------------------------------------------------------------- - SUBROUTINE setAuxiliaryVariables(num_aux,mbc,mx,my,mz,xlower,ylower,zlower,& - dxc,dyc,dzc,xyz0,xyzN,mapType,aux) - - IMPLICIT NONE - - ! Input - INTEGER, INTENT(IN) :: num_aux,mbc,mx,my,mz - REAL(kind=8), INTENT(IN) :: xlower,ylower,zlower,dxc,dyc,dzc - REAL(kind=8), DIMENSION(3), INTENT(IN) :: xyz0,xyzN - CHARACTER(20) :: mapType - - ! Output - REAL(kind=8), DIMENSION(num_aux,mx+2*mbc,my+2*mbc,mz+2*mbc), & - INTENT(OUT) :: aux - - ! Local - REAL(kind=8), DIMENSION(3) :: norm - REAL(kind=8), DIMENSION(3,4) :: poly - REAL(kind=8), DIMENSION(8) :: xcCorner,ycCorner,zcCorner,& - xpCorner,ypCorner,zpCorner - INTEGER :: i,j,k - REAL(kind=8) :: vHex,area - REAL(kind=8), DIMENSION(1) :: xc,yc,zc,xp,yp,zp - - ! Compute Auxiliary Variables - DO i = 1,mx+2*mbc - DO j = 1,my+2*mbc - DO k = 1,mz+2*mbc - - ! Compute Physical Coordinates of Corners of Hexahedron - - ! i-1/2, j-1/2, k-1/2 Corner - xcCorner(1) = xlower + (i-1)*dxc - ycCorner(1) = ylower + (j-1)*dyc - zcCorner(1) = zlower + (k-1)*dzc - xc(1) = xcCorner(1) - yc(1) = ycCorner(1) - zc(1) = zcCorner(1) - CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - xpCorner(1) = xp(1) - ypCorner(1) = yp(1) - zpCorner(1) = zp(1) - - ! i+1/2, j-1/2, k-1/2 Corner - xcCorner(2) = xcCorner(1) + dxc - ycCorner(2) = ycCorner(1) - zcCorner(2) = zcCorner(1) - xc(1) = xcCorner(2) - yc(1) = ycCorner(2) - zc(1) = zcCorner(2) - CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - xpCorner(2) = xp(1) - ypCorner(2) = yp(1) - zpCorner(2) = zp(1) - - ! i+1/2, j+1/2, k-1/2 Corner - xcCorner(3) = xcCorner(1) + dxc - ycCorner(3) = ycCorner(1) + dyc - zcCorner(3) = zcCorner(1) - xc(1) = xcCorner(3) - yc(1) = ycCorner(3) - zc(1) = zcCorner(3) - CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - xpCorner(3) = xp(1) - ypCorner(3) = yp(1) - zpCorner(3) = zp(1) - - ! i-1/2, j+1/2, k-1/2 Corner - xcCorner(4) = xcCorner(1) - ycCorner(4) = ycCorner(1) + dyc - zcCorner(4) = zcCorner(1) - xc(1) = xcCorner(4) - yc(1) = ycCorner(4) - zc(1) = zcCorner(4) - CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - xpCorner(4) = xp(1) - ypCorner(4) = yp(1) - zpCorner(4) = zp(1) - - ! i-1/2, j-1/2, k+1/2 Corner - xcCorner(5) = xcCorner(1) - ycCorner(5) = ycCorner(1) - zcCorner(5) = zcCorner(1) + dzc - xc(1) = xcCorner(5) - yc(1) = ycCorner(5) - zc(1) = zcCorner(5) - CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - xpCorner(5) = xp(1) - ypCorner(5) = yp(1) - zpCorner(5) = zp(1) - - ! i+1/2, j-1/2, k+1/2 Corner - xcCorner(6) = xcCorner(1) + dxc - ycCorner(6) = ycCorner(1) - zcCorner(6) = zcCorner(1) + dzc - xc(1) = xcCorner(6) - yc(1) = ycCorner(6) - zc(1) = zcCorner(6) - CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - xpCorner(6) = xp(1) - ypCorner(6) = yp(1) - zpCorner(6) = zp(1) - - ! i+1/2, j+1/2, k+1/2 Corner - xcCorner(7) = xcCorner(1) + dxc - ycCorner(7) = ycCorner(1) + dyc - zcCorner(7) = zcCorner(1) + dzc - xc(1) = xcCorner(7) - yc(1) = ycCorner(7) - zc(1) = zcCorner(7) - CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - xpCorner(7) = xp(1) - ypCorner(7) = yp(1) - zpCorner(7) = zp(1) - - ! i-1/2, j+1/2, k+1/2 Corner - xcCorner(8) = xcCorner(1) - ycCorner(8) = ycCorner(1) + dyc - zcCorner(8) = zcCorner(1) + dzc - xc(1) = xcCorner(8) - yc(1) = ycCorner(8) - zc(1) = zcCorner(8) - CALL mapc2pWrapper(xc,yc,zc,mz,xyz0,xyzN,mapType,xp,yp,zp) - xpCorner(8) = xp(1) - ypCorner(8) = yp(1) - zpCorner(8) = zp(1) - - ! Compute Aux Variables (Normals/Area Ratios/Volume Ratio) - - ! \xi normal - poly(:,1) = (/xpCorner(1),ypCorner(1),zpCorner(1)/) - poly(:,2) = (/xpCorner(4),ypCorner(4),zpCorner(4)/) - poly(:,3) = (/xpCorner(8),ypCorner(8),zpCorner(8)/) - poly(:,4) = (/xpCorner(5),ypCorner(5),zpCorner(5)/) - CALL unitNormal(poly(:,1),poly(:,2),poly(:,3),norm) - aux(1,i,j,k) = norm(1) ! nx(i-1/2,j,k) - aux(2,i,j,k) = norm(2) ! ny(i-1/2,j,k) - aux(3,i,j,k) = norm(3) ! nz(i-1/2,j,k) - CALL areaPolygon(poly,area) - aux(4,i,j,k) = area/(dyc*dzc) ! gamma(i-1/2,j,k) - - ! \eta normal - poly(:,1) = (/xpCorner(1),ypCorner(1),zpCorner(1)/) - poly(:,2) = (/xpCorner(5),ypCorner(5),zpCorner(5)/) - poly(:,3) = (/xpCorner(6),ypCorner(6),zpCorner(6)/) - poly(:,4) = (/xpCorner(2),ypCorner(2),zpCorner(2)/) - CALL unitNormal(poly(:,1),poly(:,2),poly(:,3),norm) - aux(5,i,j,k) = norm(1) ! nx(i,j-1/2,k) - aux(6,i,j,k) = norm(2) ! ny(i,j-1/2,k) - aux(7,i,j,k) = norm(3) ! nz(i,j-1/2,k) - CALL areaPolygon(poly,area) - aux(8,i,j,k) = area/(dxc*dzc) ! gamma(i,j-1/2,k) - - ! \phi normal - poly(:,1) = (/xpCorner(1),ypCorner(1),zpCorner(1)/) - poly(:,2) = (/xpCorner(2),ypCorner(2),zpCorner(2)/) - poly(:,3) = (/xpCorner(3),ypCorner(3),zpCorner(3)/) - poly(:,4) = (/xpCorner(4),ypCorner(4),zpCorner(4)/) - CALL unitNormal(poly(:,1),poly(:,2),poly(:,3),norm) - aux(9,i,j,k) = norm(1) ! nx(i,j,k-1/2) - aux(10,i,j,k) = norm(2) ! ny(i,j,k-1/2) - aux(11,i,j,k) = norm(3) ! nz(i,j,k-1/2) - CALL areaPolygon(poly,area) - aux(12,i,j,k) = area/(dxc*dyc) ! gamma(i,j,k-1/2) - - ! Volume of hexahedron - vHex = volumeOfHex(xpCorner,ypCorner,zpCorner) - - ! Capacity Function (Ratio of Physical:Computational volume) - aux(13,i,j,k) = vHex/(dxc*dyc*dzc) ! kappa(i,j,k) - - ! Grid Delta in Direction of Gravity - SELECT CASE(mapType) - CASE("Identity","ZeroToOne","Rotate") - aux(14,i,j,k) = ABS(zpCorner(5)-zpCorner(1)) ! dz - aux(15,i,j,k) = (zpCorner(5) + zpCorner(1))/2.d0 ! z - CASE("Spherical") - aux(14,i,j,k) = & - ABS(SQRT(xpCorner(5)**2+ypCorner(5)**2+zpCorner(5)**2) & - - SQRT(xpCorner(1)**2+ypCorner(1)**2+zpCorner(1)**2)) ! dr - aux(15,i,j,k) = (SQRT(xpCorner(5)**2+ypCorner(5)**2+zpCorner(5)**2) & - + SQRT(xpCorner(1)**2+ypCorner(1)**2+zpCorner(1)**2))/2.d0 & - - rEarth - END SELECT - - END DO - END DO - END DO - - RETURN - END SUBROUTINE setAuxiliaryVariables - -END MODULE euler3d_mappedgrid diff --git a/examples/euler_gravity_3d/rpn3_euler_mapgrid.f90 b/examples/euler_gravity_3d/rpn3_euler_mapgrid.f90 deleted file mode 100644 index 026b16998..000000000 --- a/examples/euler_gravity_3d/rpn3_euler_mapgrid.f90 +++ /dev/null @@ -1,312 +0,0 @@ -!----------------------------------------------------------------------------------- -! Name: rpn3(ixyz,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq) -! -! Description: Roe f-wave solver for the 3D Euler equations in general geometries -! -! Waves: 3 -! Equations: 5 -! -! Conserved Quantities: -! 1 density -! 2 x_momentum -! 3 y_momentum -! 4 z_momentum -! 5 energy -! -! Auxilliary Quantities: -! 1 nx xface (i-1/2,j,k) -! 2 ny xface (i-1/2,j,k) -! 3 nz xface (i-1/2,j,k) -! 4 gamma_x (area/(dyc*dzc)) -! 5 nx yface (i,j-1/2,k) -! 6 ny yface (i,j-1/2,k) -! 7 nz yface (i,j-1/2,k) -! 8 gamma_x (area/(dxc*dzc)) -! 9 nx zface (i,j,k-1/2) -! 10 ny zface (i,j,k-1/2) -! 11 nz zface (i,j,k-1/2) -! 12 gamma_z (area/(dxc*dyc)) -! 13 kappa(i,j,k) volumeOfHex/(dxc*dyc*dzc) -! 14 dz (Cartesian), dr (Spherical) -! 15 z (Cartesian), r (Spherical) -! -! Inputs: -! ixyz : direction to take slice x-direction if ixyz=1 -! y-direction if ixyz=2. -! z-direction if ixyz=3. -! maxm : max number of grid cells (less the ghost cells) -! meqn : number of equations in the system (=5) -! mwaves : nuber of waves in the system (=3) -! maux : number of auxilary variables (=15) -! mbc : number of ghost cells on either end -! mx : number of elements -! ql : state vector at left edge of each cell -! Note the i'th Riemann problem has left state qr(:,i-1) -! qr : state vector at right edge of each cell -! Note the i'th Riemann problem has right state ql(:,i) -! auxl : state of auxilary variable on left egde of cell -! auxr : state of auxilary variable on right egde of cell -! -! Outputs: -! fwave : f-wave vectors (eigenvectors of Roe matrix) -! s : eigenvalues (eigenvalues of Roe matrix) -! amdq : left-going fluctuations -! apdq : right-going fluctuations -! -! Adapted from rpn3_euler.f90 in $CLAWHOME/riemann/src -!----------------------------------------------------------------------------------- -SUBROUTINE rpn3(ixyz,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq) - - IMPLICIT NONE - - ! Input - INTEGER, INTENT(IN) :: ixyz,maxm,meqn,mwaves,maux,mbc,mx - REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(IN) :: ql,qr - REAL(kind=8), DIMENSION(maux,1-mbc:maxm+mbc), INTENT(IN) :: auxl,auxr - - ! Output - REAL(kind=8), DIMENSION(meqn,mwaves,1-mbc:maxm+mbc), INTENT(INOUT) :: fwave - REAL(kind=8), DIMENSION(mwaves,1-mbc:maxm+mbc), INTENT(INOUT) :: s - REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(INOUT) :: amdq,apdq - - ! Local Storage - REAL(kind=8), DIMENSION(meqn) :: delta,beta - REAL(kind=8), DIMENSION(-1:maxm+mbc) :: u,v,w,enth,ek,a,a2g1,psi - REAL(kind=8) :: gamma1,asqrd,pl,pr,vnl,vnr,rhsq2,rhsqrtl,rhsqrtr,nx,ny,nz,vn,& - areaRatio,dr,finterp,qinterp,quinterp,qvinterp,qwinterp - INTEGER :: i,j,mu=2,mv=3,mw=4,mws,nu,nv,nw,iratio=4,gFlux=0 - CHARACTER(20) :: nonsingularNormal - - ! Common block storage for ideal gas constant and gravity term - REAL(kind=8) :: gamma ! Ideal gas constant - REAL(kind=8) :: g_r ! Gravitational Constant Magnitude - LOGICAL :: gravity,gravityflux ! Turn Gravity Source Term On/Off - COMMON /cparam/ gamma,g_r,gravity,gravityflux - - ! Set (gamma-1) - gamma1 = gamma - 1.d0 - - ! Gravity in Energy Flux On/Off - gFlux = 0 - IF(gravityflux)gFlux=1 - - ! Make sure the number of waves is 3 for this solver - IF (mwaves /= 3) THEN - WRITE(6,*) '*** Must set mwaves=3 for this Riemann solver' - STOP - ENDIF - - ! Set mu to point to the component of the system that corresponds to momentum in - ! the direction of this slice, mv and mw to the orthogonal momentums. Set the - ! nu,nv,and nw indices for the normal vectors in a similar way. The iratio is - ! the aux array index for the gamma area scaling for mapped grids. - IF(ixyz == 1)THEN - mu = 2 - mv = 3 - mw = 4 - iratio = 4 - ELSE IF(ixyz == 2)THEN - mu = 3 - mv = 4 - mw = 2 - iratio = 8 - ELSE IF (ixyz == 3)THEN - mu = 4 - mv = 2 - mw = 3 - iratio = 12 - ELSE - WRITE(*,*)"Something is wrong with ixyz..." - STOP - ENDIF - nu = (ixyz-1)*4 + mu - 1 - nv = (ixyz-1)*4 + mv - 1 - nw = (ixyz-1)*4 + mw - 1 - - ! Note that notation for u,v, and w reflects assumption that the - ! Riemann problems are in the x-direction with u in the normal - ! direction and v and w in the orthogonal directions, but with the - ! above definitions of mu, mv, and mw the routine also works with - ! ixyz=2 and ixyz = 3 - ! and returns, for example, f0 as the Godunov flux g0 for the - ! Riemann problems u_t + g(u)_y = 0 in the y-direction. - - ! Initialize fwaves and speeds - fwave = 0.d0 - s = 0.d0 - - ! Loop over grid cell interfaces - DO i = 2-mbc, mx+mbc - IF (qr(1,i-1) <= 0.d0 .OR. ql(1,i) <= 0.d0) THEN - WRITE(*,*) i, mu, mv, mw - WRITE(*,'(5e12.4)') (qr(j,i-1),j=1,5) - WRITE(*,'(5e12.4)') (ql(j,i),j=1,5) - IF (ixyz == 1) WRITE(6,*) '*** (rpn3) rho <= 0 in x-sweep at ',i - IF (ixyz == 2) WRITE(6,*) '*** (rpn3) rho <= 0 in y-sweep at ',i - IF (ixyz == 3) WRITE(6,*) '*** (rpn3) rho <= 0 in z-sweep at ',i - WRITE(6,*) 'stopped with rho <= 0...' - STOP - ENDIF - rhsqrtl = SQRT(qr(1,i-1)) - rhsqrtr = SQRT(ql(1,i)) - pl = gamma1*(qr(5,i-1) - 0.5d0*(qr(mu,i-1)**2 + & - qr(mv,i-1)**2 + qr(mw,i-1)**2)/qr(1,i-1) & - - qr(1,i-1)*(g_r)*auxr(15,i-1)*gFlux ) - pr = gamma1*(ql(5,i) - 0.5d0*(ql(mu,i)**2 + & - ql(mv,i)**2 + ql(mw,i)**2)/ql(1,i) & - - ql(1,i)*(g_r)*auxl(15,i)*gFlux ) - rhsq2 = rhsqrtl + rhsqrtr - u(i) = (qr(mu,i-1)/rhsqrtl + ql(mu,i)/rhsqrtr) / rhsq2 - v(i) = (qr(mv,i-1)/rhsqrtl + ql(mv,i)/rhsqrtr) / rhsq2 - w(i) = (qr(mw,i-1)/rhsqrtl + ql(mw,i)/rhsqrtr) / rhsq2 - enth(i) = ((qr(5,i-1)+pl)/rhsqrtl + (ql(5,i)+pr)/rhsqrtr) / rhsq2 - psi(i) = (qr(1,i-1)*(g_r)*auxr(15,i-1)/rhsqrtl & - + ql(1,i)*(g_r)*auxl(15,i)/rhsqrtr)/rhsq2*gFlux - ek(i) = 0.5d0*(u(i)**2 + v(i)**2 + w(i)**2) - asqrd = gamma1*(enth(i) - ek(i)) - IF (i>=0 .AND. i<=mx .AND. asqrd <= 0.d0) THEN - IF (ixyz == 1) WRITE(6,*) '*** (rpn3) a**2 <= 0 in x-sweep at ',i - IF (ixyz == 2) WRITE(6,*) '*** (rpn3) a**2 <= 0 in y-sweep at ',i - IF (ixyz == 3) WRITE(6,*) '*** (rpn3) a**2 <= 0 in z-sweep at ',i - WRITE(6,*) 'stopped with a**2 < 0...' - STOP - ENDIF - a(i) = SQRT(asqrd) - a2g1(i) = asqrd/gamma1 - END DO - - ! Now split the jump in q1d at each interface into f-waves - ! find alpha(1) thru alpha(5), the coefficients of the 5 eigenvectors: - DO i = 2-mbc, mx+mbc - - ! Normal Vectors, Normal Velocity, and Area Ratio - nx = auxl(nu,i) - ny = auxl(nv,i) - nz = auxl(nw,i) - vn = u(i)*nx + v(i)*ny + w(i)*nz - areaRatio = auxl(iratio,i) - - ! Compute the jumps (delta F) in the flux functions - pl = gamma1*(ql(5,i) & - - 0.5d0*(ql(mu,i)**2+ql(mv,i)**2+ql(mw,i)**2)/ql(1,i) & - - ql(1,i)*(g_r)*auxl(15,i)*gFlux) - pr = gamma1*(qr(5,i-1) & - - 0.5d0*(qr(mu,i-1)**2+qr(mv,i-1)**2+qr(mw,i-1)**2)/qr(1,i-1) & - - qr(1,i-1)*(g_r)*auxr(15,i-1)*gFlux) - vnl = (ql(mu,i)*nx + ql(mv,i)*ny + ql(mw,i)*nz)/ql(1,i) - vnr = (qr(mu,i-1)*nx + qr(mv,i-1)*ny + qr(mw,i-1)*nz)/qr(1,i-1) - delta(1) = ql(1,i)*vnl - qr(1,i-1)*vnr - delta(2) = (vnl*ql(mu,i) + pl*nx) - (vnr*qr(mu,i-1) + pr*nx) - delta(3) = (vnl*ql(mv,i) + pl*ny) - (vnr*qr(mv,i-1) + pr*ny) - delta(4) = (vnl*ql(mw,i) + pl*nz) - (vnr*qr(mw,i-1) + pr*nz) - delta(5) = ((ql(5,i) + pl)*vnl) - ((qr(5,i-1) + pr)*vnr) - - ! Modify delta(2/3/4/5) to account for gravitational term - IF(gravity .AND. ixyz==3)THEN - dr = (auxl(14,i) + auxl(14,i-1))*0.5d0 - finterp = auxl(14,i-1)/(auxl(14,i)+auxl(14,i-1)) - qinterp = ql(1,i)*finterp + qr(1,i-1)*(1.d0-finterp) - quinterp = ql(mu,i)*finterp + qr(mu,i-1)*(1.d0-finterp) - qvinterp = ql(mv,i)*finterp + qr(mv,i-1)*(1.d0-finterp) - qwinterp = ql(mw,i)*finterp + qr(mw,i-1)*(1.d0-finterp) - delta(2) = delta(2) + dr*g_r*nx*qinterp - delta(3) = delta(3) + dr*g_r*ny*qinterp - delta(4) = delta(4) + dr*g_r*nz*qinterp - IF(.NOT.gravityflux)& - delta(5) = delta(5) + dr*g_r*(quinterp*nx + qvinterp*ny + qwinterp*nz) - END IF - - ! Pick the largest normal to avoid using a singular direction - IF(ABS(nx) >= ABS(ny) .AND. ABS(nx) >= ABS(nz))THEN - nonsingularNormal = "x" - ELSEIF(ABS(ny) >= ABS(nx) .AND. ABS(ny) >= ABS(nz))THEN - nonsingularNormal = "y" - ELSEIF(ABS(nz) >= ABS(nx) .AND. ABS(nz) >= ABS(ny))THEN - nonsingularNormal = "z" - ELSE - WRITE(*,*)"Invalid Normal Vector..." - WRITE(*,*)"nx: ",nx,"ny: ",ny,"nz: ",nz - STOP - END IF - - beta(4) = ((enth(i) + psi(i) - 2.d0*ek(i))*delta(1) & - + u(i)*delta(2) + v(i)*delta(3) + w(i)*delta(4) - delta(5))/a2g1(i) - beta(5) = ( (a(i)-vn)*delta(1) - a(i)*beta(4) & - + nx*delta(2) + ny*delta(3) + nz*delta(4) )/(2.d0*a(i)) - beta(1) = delta(1) - beta(4) - beta(5) - - SELECT CASE(TRIM(nonsingularNormal)) - CASE("x") - beta(2) = ( (v(i)-vn*ny)*delta(1) + nx*ny*delta(2) & - + (ny**2-1.d0)*delta(3) + ny*nz*delta(4) )/nx - beta(3) = (-(w(i)-vn*nz)*delta(1) - nx*nz*delta(2) & - - ny*nz*delta(3) - (nz**2-1.d0)*delta(4) )/nx - CASE("y") - beta(2) = (-(u(i)-vn*nx)*delta(1) - (nx**2-1.d0)*delta(2) & - - ny*nx*delta(3) - nx*nz*delta(4) )/ny - beta(3) = ( (w(i)-vn*nz)*delta(1) + nx*nz*delta(2) & - + ny*nz*delta(3) + (nz**2-1.d0)*delta(4) )/ny - CASE("z") - beta(2) = ( (u(i)-vn*nx)*delta(1) + (nx**2-1.d0)*delta(2) & - + nx*ny*delta(3) + nz*nx*delta(4) )/nz - beta(3) = (-(v(i)-vn*ny)*delta(1) - nx*ny*delta(2) & - - (ny**2-1.d0)*delta(3) - nz*ny*delta(4) )/nz - END SELECT - - ! Compute the f-waves. - ! Note that the 2-wave, 3-wave and 4-wave travel at the same speed - ! and are lumped together in fwave(.,2,.). The 5-wave is then stored - ! in fwave(.,3,.). - fwave(1,1,i) = beta(1) - fwave(mu,1,i) = beta(1)*(u(i)-a(i)*nx) - fwave(mv,1,i) = beta(1)*(v(i)-a(i)*ny) - fwave(mw,1,i) = beta(1)*(w(i)-a(i)*nz) - fwave(5,1,i) = beta(1)*(enth(i) + psi(i) - a(i)*vn) - s(1,i) = (vn-a(i))*areaRatio - - fwave(1,2,i) = beta(4) - SELECT CASE(TRIM(nonsingularNormal)) - CASE("x") - fwave(mu,2,i) = beta(4)*u(i) + beta(2)*ny - beta(3)*nz - fwave(mv,2,i) = beta(4)*v(i) - beta(2)*nx - fwave(mw,2,i) = beta(4)*w(i) + beta(3)*nx - fwave(5,2,i) = beta(4)*(ek(i) + psi(i)) & - + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(w(i)*nx-u(i)*nz) - CASE("y") - fwave(mu,2,i) = beta(4)*u(i) + beta(2)*ny - fwave(mv,2,i) = beta(4)*v(i) - beta(2)*nx + beta(3)*nz - fwave(mw,2,i) = beta(4)*w(i) - beta(3)*ny - fwave(5,2,i) = beta(4)*(ek(i) + psi(i)) & - + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(v(i)*nz-w(i)*ny) - CASE("z") - fwave(mu,2,i) = beta(4)*u(i) - beta(2)*nz - fwave(mv,2,i) = beta(4)*v(i) + beta(3)*nz - fwave(mw,2,i) = beta(4)*w(i) + beta(2)*nx - beta(3)*ny - fwave(5,2,i) = beta(4)*(ek(i) + psi(i)) & - + beta(2)*(w(i)*nx-u(i)*nz) + beta(3)*(v(i)*nz-w(i)*ny) - END SELECT - s(2,i) = vn*areaRatio - - fwave(1,3,i) = beta(5) - fwave(mu,3,i) = beta(5)*(u(i)+a(i)*nx) - fwave(mv,3,i) = beta(5)*(v(i)+a(i)*ny) - fwave(mw,3,i) = beta(5)*(w(i)+a(i)*nz) - fwave(5,3,i) = beta(5)*(enth(i) + psi(i) + a(i)*vn) - s(3,i) = (vn+a(i))*areaRatio - - ! Compute fluctuations amdq and apdq - ! amdq = SUM s*fwave over left-going waves - ! apdq = SUM s*fwave over right-going waves - amdq(:,i) = 0.d0 - apdq(:,i) = 0.d0 - DO mws=1,mwaves - fwave(:,mws,i) = fwave(:,mws,i)*areaRatio - IF (s(mws,i) < 0.d0) THEN - amdq(:,i) = amdq(:,i) + fwave(:,mws,i) - ELSE - apdq(:,i) = apdq(:,i) + fwave(:,mws,i) - ENDIF - END DO - - END DO -END SUBROUTINE rpn3 diff --git a/examples/euler_gravity_3d/rpt3_euler_mapgrid.f90 b/examples/euler_gravity_3d/rpt3_euler_mapgrid.f90 deleted file mode 100644 index 09f215eaf..000000000 --- a/examples/euler_gravity_3d/rpt3_euler_mapgrid.f90 +++ /dev/null @@ -1,295 +0,0 @@ -!----------------------------------------------------------------------------------- -! Name: rpt3(ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq, -! bmasdq,bpasdq) -! -! Description: Roe f-wave solver in the transverse direction for the -! 3D Euler equations in general geometries -! -! Inputs: -! ixyz : direction to take slice x-direction if ixyz=1 -! y-direction if ixyz=2. -! z-direction if ixyz=3. -! icoor : direction in which the transverse solve should be -! performed. icoor=2 (split in y-like direction) -! icoor=3 (split in z-like direction) -! ixyz=1, icoor=2 means asdq=A^*\Dq, split in y into: -! bmasdq = B^-A^*\Dq, -! bpasdq = B^+A^*\Dq. -! ixyz=2, icoor=2 means asdq=B^*\Dq, split in z into: -! bmasdq = C^-B^*\Dq, -! bpasdq = C^+B^*\Dq. -! imp : index for aux arrays -! maxm : max number of grid cells (less the ghost cells) -! meqn : number of equations in the system -! mwaves : nuber of waves in the system -! maux : number of auxilary equations -! mbc : number of ghost cells on either end -! mx : number of elements -! ql : state vector at left edge of each cell -! Note the i'th Riemann problem has left state qr(:,i-1) -! qr : state vector at right edge of each cell -! Note the i'th Riemann problem has right state ql(:,i) -! aux1 : -! aux2 : -! aux3 : -! asdq : array of flux differences (A*\Dq). asdq(:,i) is the -! flux difference propagating away from the interface -! between cells i-1 and i. Note: asdq represents -! B*\Dq if ixyz=2 or C*\Dq if ixyz=3. -! -! Outputs: -! bmasdq : left-going flux differences -! bpasdq : right-going flux differences -! -! Adapted from rpt3_euler.f90 in $CLAWHOME/riemann/src -!----------------------------------------------------------------------------------- -SUBROUTINE rpt3(ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,& - asdq,bmasdq,bpasdq) - - IMPLICIT NONE - - ! Input - INTEGER, INTENT(IN) :: ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx - REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(IN) :: ql,qr,asdq - REAL(kind=8), DIMENSION(maux,1-mbc:maxm+mbc,3), INTENT(IN) :: aux1,aux2,aux3 - - ! Output - REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(INOUT) :: bmasdq,bpasdq - - ! Local Storage - REAL(kind=8), DIMENSION(meqn) :: beta - REAL(kind=8), DIMENSION(-1:maxm+mbc) :: u2v2w2,u,v,w,enth,a,g1a2,euv - REAL(kind=8) :: gamma1,asqrd,pl,pr,rhsq2,rhsqrtl,rhsqrtr,nx,ny,nz,vn,& - areaRatioMinus,areaRatioPlus - REAL(kind=8), DIMENSION(meqn,mwaves) :: waveb - REAL(kind=8), DIMENSION(mwaves) :: sb - INTEGER :: i,j,mu=2,mv=3,mw=4,mws,nu=1,nv=2,nw=3,iratio=4 - CHARACTER(20) :: nonsingularNormal - - ! Common block storage for ideal gas constant - REAL(kind=8) :: gamma - COMMON /cparam/ gamma - - ! Set (gamma-1) - gamma1 = gamma - 1.d0 - - ! Make sure the number of waves is 3 for the 3D Euler equations - IF (mwaves /= 3) THEN - WRITE(6,*) '*** Should have mwaves=3 for this Riemann solver' - STOP - ENDIF - - ! Set mu to point to the component of the system that corresponds to momentum in - ! the direction of this slice, mv and mw to the orthogonal momentums. - ! - ! ixyz indicates the direction of the original Riemann solve, - ! called the x-like direction in the table below: - ! - ! x-like direction y-like direction z-like direction - ! ixyz=1: x y z - ! ixyz=2: y z x - ! ixyz=3: z x y - IF(ixyz == 1)THEN ! x-like direction - mu = 2 - mv = 3 - mw = 4 - ! Choose the coordinate direction to solve the Riemann problem - ! Permute nu,nv,nw for the original ixyz incident direction (x-like) - IF(icoor == 2)THEN ! y-like direction - nu = 5 - nv = 6 - nw = 7 - iratio = 8 - ELSEIF(icoor == 3)THEN ! z-like direction - nu = 9 - nv = 10 - nw = 11 - iratio = 12 - END IF - ELSEIF(ixyz == 2)THEN ! y-like direction - mu = 3 - mv = 4 - mw = 2 - ! Choose the coordinate direction to solve the Riemann problem - ! Permute nu,nv,nw for the original ixyz incident direction (y-like) - IF(icoor == 2)THEN ! z-like direction - nu = 10 - nv = 11 - nw = 9 - iratio = 12 - ELSEIF(icoor == 3)THEN ! x-like direction - nu = 2 - nv = 3 - nw = 1 - iratio = 4 - END IF - ELSEIF(ixyz == 3)THEN ! z-like direction - mu = 4 - mv = 2 - mw = 3 - ! Choose the coordinate direction to solve the Riemann problem - ! Permute nu,nv,nw for the original ixyz incident direction (z-like) - IF(icoor == 2)THEN ! x-like direction - nu = 3 - nv = 1 - nw = 2 - iratio = 4 - ELSEIF(icoor == 3)THEN ! y-like direction - nu = 7 - nv = 5 - nw = 6 - iratio = 8 - END IF - ENDIF - - ! Note that notation for u,v, and w reflects assumption that the - ! Riemann problems are in the x-direction with u in the normal - ! direction and v and w in the orthogonal directions, but with the - ! above definitions of mu, mv, and mw the routine also works with - ! ixyz=2 and ixyz = 3 - ! and returns, for example, f0 as the Godunov flux g0 for the - ! Riemann problems u_t + g(u)_y = 0 in the y-direction. - ! - ! Compute the Roe-averaged variables needed in the Roe solver. - - ! Loop over grid cell interfaces - DO i = 2-mbc, mx+mbc - IF (qr(1,i-1) <= 0.d0 .OR. ql(1,i) <= 0.d0) THEN - WRITE(*,*) i, mu, mv, mw - WRITE(*,'(5e12.4)') (qr(j,i-1),j=1,5) - WRITE(*,'(5e12.4)') (ql(j,i),j=1,5) - IF (ixyz == 1) WRITE(6,*) '*** (rpt3) rho <= 0 in x-sweep at ',i - IF (ixyz == 2) WRITE(6,*) '*** (rpt3) rho <= 0 in y-sweep at ',i - IF (ixyz == 3) WRITE(6,*) '*** (rpt3) rho <= 0 in z-sweep at ',i - WRITE(6,*) 'stopped with rho <= 0...' - STOP - ENDIF - rhsqrtl = SQRT(qr(1,i-1)) - rhsqrtr = SQRT(ql(1,i)) - pl = gamma1*(qr(5,i-1) - 0.5d0*(qr(mu,i-1)**2 + & - qr(mv,i-1)**2 + qr(mw,i-1)**2)/qr(1,i-1)) - pr = gamma1*(ql(5,i) - 0.5d0*(ql(mu,i)**2 + & - ql(mv,i)**2 + ql(mw,i)**2)/ql(1,i)) - rhsq2 = rhsqrtl + rhsqrtr - u(i) = (qr(mu,i-1)/rhsqrtl + ql(mu,i)/rhsqrtr) / rhsq2 - v(i) = (qr(mv,i-1)/rhsqrtl + ql(mv,i)/rhsqrtr) / rhsq2 - w(i) = (qr(mw,i-1)/rhsqrtl + ql(mw,i)/rhsqrtr) / rhsq2 - enth(i) = (((qr(5,i-1)+pl)/rhsqrtl & - + (ql(5,i)+pr)/rhsqrtr)) / rhsq2 - u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2 - asqrd = gamma1*(enth(i) - .5d0*u2v2w2(i)) - - IF (i>=0 .AND. i<=mx .AND. asqrd <= 0.d0) THEN - IF (ixyz == 1) WRITE(6,*) '*** (rpt3) a**2 <= 0 in x-sweep at ',i - IF (ixyz == 2) WRITE(6,*) '*** (rpt3) a**2 <= 0 in y-sweep at ',i - IF (ixyz == 3) WRITE(6,*) '*** (rpt3) a**2 <= 0 in z-sweep at ',i - WRITE(6,*) 'stopped with a**2 < 0...' - STOP - ENDIF - a(i) = SQRT(asqrd) - g1a2(i) = gamma1 / asqrd - euv(i) = enth(i) - u2v2w2(i) - END DO - - ! Solve the Riemann Problem - DO i = 2-mbc, mx+mbc - - nx = aux2(nu,i,2) - ny = aux2(nv,i,2) - nz = aux2(nw,i,2) - vn = u(i)*nx + v(i)*ny + w(i)*nz - IF(icoor == 2)THEN - areaRatioMinus = aux2(iratio,i-2+imp,2) - areaRatioPlus = aux3(iratio,i-2+imp,2) - ELSEIF(icoor == 3)THEN - areaRatioMinus = aux2(iratio,i-2+imp,2) - areaRatioPlus = aux2(iratio,i-2+imp,3) - ENDIF - ! Pick the largest normal to avoid using a singular direction - IF(ABS(nx) >= ABS(ny) .AND. ABS(nx) >= ABS(nz))THEN - nonsingularNormal = "x" - ELSEIF(ABS(ny) >= ABS(nx) .AND. ABS(ny) >= ABS(nz))THEN - nonsingularNormal = "y" - ELSEIF(ABS(nz) >= ABS(nx) .AND. ABS(nz) >= ABS(ny))THEN - nonsingularNormal = "z" - ELSE - WRITE(*,*)"Invalid Normal Vector..." - WRITE(*,*)"nx: ",nx,"ny: ",ny,"nz: ",nz - STOP - END IF - - beta(4) = g1a2(i)*(euv(i)*asdq(1,i) & - + u(i)*asdq(mu,i) + v(i)*asdq(mv,i) + w(i)*asdq(mw,i) - asdq(5,i)) - beta(5) = ( (a(i)-vn)*asdq(1,i) - a(i)*beta(4) & - + nx*asdq(mu,i) + ny*asdq(mv,i) + nz*asdq(mw,i) )/(2.d0*a(i)) - beta(1) = asdq(1,i) - beta(4) - beta(5) - - SELECT CASE(TRIM(nonsingularNormal)) - CASE("x") - beta(2) = ( (v(i)-vn*ny)*asdq(1,i) + nx*ny*asdq(mu,i) & - + (ny**2-1)*asdq(mv,i) + ny*nz*asdq(mw,i) )/nx - beta(3) = (-(w(i)-vn*nz)*asdq(1,i) - nx*nz*asdq(mu,i) & - - ny*nz*asdq(mv,i) - (nz**2-1)*asdq(mw,i) )/nx - CASE("y") - beta(2) = (-(u(i)-vn*nx)*asdq(1,i) - (nx**2-1)*asdq(mu,i) & - - ny*nx*asdq(mv,i) - nx*nz*asdq(mw,i) )/ny - beta(3) = ( (w(i)-vn*nz)*asdq(1,i) + nx*nz*asdq(mu,i) & - + ny*nz*asdq(mv,i) + (nz**2-1)*asdq(mw,i) )/ny - CASE("z") - beta(2) = ( (u(i)-vn*nx)*asdq(1,i) + (nx**2-1)*asdq(mu,i) & - + nx*ny*asdq(mv,i) + nz*nx*asdq(mw,i) )/nz - beta(3) = (-(v(i)-vn*ny)*asdq(1,i) - nx*ny*asdq(mu,i) & - - (ny**2-1)*asdq(mv,i) - nz*ny*asdq(mw,i) )/nz - END SELECT - - waveb(1,1) = beta(1) - waveb(mu,1) = beta(1)*(u(i)-a(i)*nx) - waveb(mv,1) = beta(1)*(v(i)-a(i)*ny) - waveb(mw,1) = beta(1)*(w(i)-a(i)*nz) - waveb(5,1) = beta(1)*(enth(i) - vn*a(i)) - sb(1) = vn - a(i) - - waveb(1,2) = beta(4) - SELECT CASE(TRIM(nonsingularNormal)) - CASE("x") - waveb(mu,2) = beta(4)*u(i) + beta(2)*ny - beta(3)*nz - waveb(mv,2) = beta(4)*v(i) - beta(2)*nx - waveb(mw,2) = beta(4)*w(i) + beta(3)*nx - waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & - + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(w(i)*nx-u(i)*nz) - CASE("y") - waveb(mu,2) = beta(4)*u(i) + beta(2)*ny - waveb(mv,2) = beta(4)*v(i) - beta(2)*nx + beta(3)*nz - waveb(mw,2) = beta(4)*w(i) - beta(3)*ny - waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & - + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(v(i)*nz-w(i)*ny) - CASE("z") - waveb(mu,2) = beta(4)*u(i) - beta(2)*nz - waveb(mv,2) = beta(4)*v(i) + beta(3)*nz - waveb(mw,2) = beta(4)*w(i) + beta(2)*nx - beta(3)*ny - waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & - + beta(2)*(w(i)*nx-u(i)*nz) + beta(3)*(v(i)*nz-w(i)*ny) - END SELECT - sb(2) = vn - - waveb(1,3) = beta(5) - waveb(mu,3) = beta(5)*(u(i) + a(i)*nx) - waveb(mv,3) = beta(5)*(v(i) + a(i)*ny) - waveb(mw,3) = beta(5)*(w(i) + a(i)*nz) - waveb(5,3) = beta(5)*(enth(i) + vn*a(i)) - sb(3) = vn + a(i) - - bmasdq(:,i) = 0.d0 - bpasdq(:,i) = 0.d0 - DO mws = 1,mwaves - IF (sb(mws) < 0.d0)THEN - waveb(:,mws) = waveb(:,mws)*areaRatioMinus - bmasdq(:,i) = bmasdq(:,i) + sb(mws)*waveb(:,mws) - ELSE - waveb(:,mws) = waveb(:,mws)*areaRatioPlus - bpasdq(:,i) = bpasdq(:,i) + sb(mws)*waveb(:,mws) - END IF - END DO - END DO - -END SUBROUTINE rpt3 diff --git a/examples/euler_gravity_3d/rptt3_euler_mapgrid.f90 b/examples/euler_gravity_3d/rptt3_euler_mapgrid.f90 deleted file mode 100644 index 9a4392769..000000000 --- a/examples/euler_gravity_3d/rptt3_euler_mapgrid.f90 +++ /dev/null @@ -1,306 +0,0 @@ -!----------------------------------------------------------------------------------- -! Name: rptt3(ixyz,icoor,imp,impt,maxm,meqn,mwaves,maux,mbc,mx,ql,qr, -! aux1,aux2,aux3,asdq,cmbsasdq,cpbsasdq) -! -! Description: Roe f-wave solver in the transverse direction for the -! 3D Euler equations in general geometries -! -! Inputs: -! ixyz : direction to take slice x-direction if ixyz=1 -! y-direction if ixyz=2. -! z-direction if ixyz=3. -! icoor : direction in which the transverse solve should be -! performed. icoor=2 (split in y-like direction) -! icoor=3 (split in z-like direction) -! ixyz=1, icoor=3 means bsasdq=B*A*\Dq, split in z into: -! cmbsasdq = C^-B*A*\Dq, -! cpbsasdq = C^+B*A*\Dq. -! ixyz=2, icoor=3 means bsasdq=C*B*\Dq, split in x into: -! cmbsasdq = A^-C*B*\Dq, -! cpbsasdq = A^+C*B*\Dq. -! imp : index for aux arrays -! impt : index for aux arrays -! maxm : max number of grid cells (less the ghost cells) -! meqn : number of equations in the system -! mwaves : nuber of waves in the system -! maux : number of auxilary equations -! mbc : number of ghost cells on either end -! mx : number of elements -! ql : state vector at left edge of each cell -! Note the i'th Riemann problem has left state qr(:,i-1) -! qr : state vector at right edge of each cell -! Note the i'th Riemann problem has right state ql(:,i) -! aux1 : -! aux2 : -! aux3 : -! asdq : array of flux differences (A*\Dq). asdq(:,i) is the -! flux difference propagating away from the interface -! between cells i-1 and i. Note: asdq represents -! B*\Dq if ixyz=2 or C*\Dq if ixyz=3. -! -! Outputs: -! cmbsasdq : left-going flux differences -! cpbsasdq : right-going flux differences -! -! Adapted from rptt3_euler.f90 in $CLAWHOME/riemann/src -!----------------------------------------------------------------------------------- -SUBROUTINE rptt3(ixyz,icoor,imp,impt,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,& - aux1,aux2,aux3,bsasdq,cmbsasdq,cpbsasdq) - - IMPLICIT NONE - - ! Input - INTEGER, INTENT(IN) :: ixyz,icoor,imp,impt,maxm,meqn,mwaves,maux,mbc,mx - REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(IN) :: ql,qr,bsasdq - REAL(kind=8), DIMENSION(maux,1-mbc:maxm+mbc,3), INTENT(IN) :: aux1,aux2,aux3 - - ! Output - REAL(kind=8), DIMENSION(meqn,1-mbc:maxm+mbc), INTENT(INOUT) :: cmbsasdq,cpbsasdq - - ! Local Storage - REAL(kind=8), DIMENSION(meqn) :: beta - REAL(kind=8), DIMENSION(-1:maxm+mbc) :: u2v2w2,u,v,w,enth,a,g1a2,euv - REAL(kind=8) :: gamma1,asqrd,pl,pr,rhsq2,rhsqrtl,rhsqrtr,nx,ny,nz,vn,& - areaRatioMinus,areaRatioPlus - REAL(kind=8), DIMENSION(meqn,mwaves) :: waveb - REAL(kind=8), DIMENSION(mwaves) :: sb - INTEGER :: i,j,mu=2,mv=3,mw=4,mws,nu=1,nv=2,nw=3,iratio=4 - CHARACTER(20) :: nonsingularNormal - - ! Common block storage for ideal gas constant - REAL(kind=8) :: gamma - COMMON /cparam/ gamma - - ! Set (gamma-1) - gamma1 = gamma - 1.d0 - - ! Make sure the number of waves is 3 for the 3D Euler equations - IF (mwaves /= 3) THEN - WRITE(6,*) '*** Should have mwaves=3 for this Riemann solver' - STOP - ENDIF - - ! Set mu to point to the component of the system that corresponds to momentum in - ! the direction of this slice, mv and mw to the orthogonal momentums. - ! - ! ixyz indicates the direction of the original Riemann solve, - ! called the x-like direction in the table below: - ! - ! x-like direction y-like direction z-like direction - ! ixyz=1: x y z - ! ixyz=2: y z x - ! ixyz=3: z x y - IF(ixyz == 1)THEN ! x-like direction - mu = 2 - mv = 3 - mw = 4 - ! Choose the coordinate direction to solve the Riemann problem - ! Permute nu,nv,nw for the original ixyz incident direction (x-like) - IF(icoor == 2)THEN ! y-like direction - nu = 5 - nv = 6 - nw = 7 - iratio = 8 - ELSEIF(icoor == 3)THEN ! z-like direction - nu = 9 - nv = 10 - nw = 11 - iratio = 12 - END IF - ELSEIF(ixyz == 2)THEN ! y-like direction - mu = 3 - mv = 4 - mw = 2 - ! Choose the coordinate direction to solve the Riemann problem - ! Permute nu,nv,nw for the original ixyz incident direction (y-like) - IF(icoor == 2)THEN ! z-like direction - nu = 10 - nv = 11 - nw = 9 - iratio = 12 - ELSEIF(icoor == 3)THEN ! x-like direction - nu = 2 - nv = 3 - nw = 1 - iratio = 4 - END IF - ELSEIF(ixyz == 3)THEN ! z-like direction - mu = 4 - mv = 2 - mw = 3 - ! Choose the coordinate direction to solve the Riemann problem - ! Permute nu,nv,nw for the original ixyz incident direction (z-like) - IF(icoor == 2)THEN ! x-like direction - nu = 3 - nv = 1 - nw = 2 - iratio = 4 - ELSEIF(icoor == 3)THEN ! y-like direction - nu = 7 - nv = 5 - nw = 6 - iratio = 8 - END IF - ENDIF - - ! Note that notation for u,v, and w reflects assumption that the - ! Riemann problems are in the x-direction with u in the normal - ! direction and v and w in the orthogonal directions, but with the - ! above definitions of mu, mv, and mw the routine also works with - ! ixyz=2 and ixyz = 3 - ! and returns, for example, f0 as the Godunov flux g0 for the - ! Riemann problems u_t + g(u)_y = 0 in the y-direction. - ! - ! Compute the Roe-averaged variables needed in the Roe solver. - - ! Loop over grid cell interfaces - DO i = 2-mbc, mx+mbc - IF (qr(1,i-1) <= 0.d0 .OR. ql(1,i) <= 0.d0) THEN - WRITE(*,*) i, mu, mv, mw - WRITE(*,'(5e12.4)') (qr(j,i-1),j=1,5) - WRITE(*,'(5e12.4)') (ql(j,i),j=1,5) - IF (ixyz == 1) WRITE(6,*) '*** (rptt3) rho <= 0 in x-sweep at ',i - IF (ixyz == 2) WRITE(6,*) '*** (rptt3) rho <= 0 in y-sweep at ',i - IF (ixyz == 3) WRITE(6,*) '*** (rptt3) rho <= 0 in z-sweep at ',i - WRITE(6,*) 'stopped with rho <= 0...' - STOP - ENDIF - rhsqrtl = SQRT(qr(1,i-1)) - rhsqrtr = SQRT(ql(1,i)) - pl = gamma1*(qr(5,i-1) - 0.5d0*(qr(mu,i-1)**2 + & - qr(mv,i-1)**2 + qr(mw,i-1)**2)/qr(1,i-1)) - pr = gamma1*(ql(5,i) - 0.5d0*(ql(mu,i)**2 + & - ql(mv,i)**2 + ql(mw,i)**2)/ql(1,i)) - rhsq2 = rhsqrtl + rhsqrtr - u(i) = (qr(mu,i-1)/rhsqrtl + ql(mu,i)/rhsqrtr) / rhsq2 - v(i) = (qr(mv,i-1)/rhsqrtl + ql(mv,i)/rhsqrtr) / rhsq2 - w(i) = (qr(mw,i-1)/rhsqrtl + ql(mw,i)/rhsqrtr) / rhsq2 - enth(i) = (((qr(5,i-1)+pl)/rhsqrtl & - + (ql(5,i)+pr)/rhsqrtr)) / rhsq2 - u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2 - asqrd = gamma1*(enth(i) - .5d0*u2v2w2(i)) - - IF (i>=0 .AND. i<=mx .AND. asqrd <= 0.d0) THEN - IF (ixyz == 1) WRITE(6,*) '*** (rptt3) a**2 <= 0 in x-sweep at ',i - IF (ixyz == 2) WRITE(6,*) '*** (rptt3) a**2 <= 0 in y-sweep at ',i - IF (ixyz == 3) WRITE(6,*) '*** (rptt3) a**2 <= 0 in z-sweep at ',i - WRITE(6,*) 'stopped with a**2 < 0...' - STOP - ENDIF - a(i) = SQRT(asqrd) - g1a2(i) = gamma1 / asqrd - euv(i) = enth(i) - u2v2w2(i) - END DO - - ! Solve the Riemann Problem - DO i = 2-mbc, mx+mbc - - nx = aux2(nu,i,2) - ny = aux2(nv,i,2) - nz = aux2(nw,i,2) - vn = u(i)*nx + v(i)*ny + w(i)*nz - IF(icoor == 2)THEN - IF(impt==1)THEN - areaRatioMinus = aux2(iratio,i-2+imp,1) - areaRatioPlus = aux3(iratio,i-2+imp,1) - ELSEIF(impt==2)THEN - areaRatioMinus = aux2(iratio,i-2+imp,3) - areaRatioPlus = aux3(iratio,i-2+imp,3) - END IF - ELSEIF(icoor == 3)THEN - IF(impt == 1)THEN - areaRatioMinus = aux1(iratio,i-2+imp,2) - areaRatioPlus = aux1(iratio,i-2+imp,3) - ELSEIF(impt == 2)THEN - areaRatioMinus = aux3(iratio,i-2+imp,2) - areaRatioPlus = aux3(iratio,i-2+imp,3) - END IF - END IF - ! Pick the largest normal to avoid using a singular direction - IF(ABS(nx) >= ABS(ny) .AND. ABS(nx) >= ABS(nz))THEN - nonsingularNormal = "x" - ELSEIF(ABS(ny) >= ABS(nx) .AND. ABS(ny) >= ABS(nz))THEN - nonsingularNormal = "y" - ELSEIF(ABS(nz) >= ABS(nx) .AND. ABS(nz) >= ABS(ny))THEN - nonsingularNormal = "z" - ELSE - WRITE(*,*)"Invalid Normal Vector..." - WRITE(*,*)"nx: ",nx,"ny: ",ny,"nz: ",nz - STOP - END IF - - beta(4) = g1a2(i)*(euv(i)*bsasdq(1,i) & - + u(i)*bsasdq(mu,i) + v(i)*bsasdq(mv,i) + w(i)*bsasdq(mw,i) - bsasdq(5,i)) - beta(5) = ( (a(i)-vn)*bsasdq(1,i) - a(i)*beta(4) & - + nx*bsasdq(mu,i) + ny*bsasdq(mv,i) + nz*bsasdq(mw,i) )/(2.d0*a(i)) - beta(1) = bsasdq(1,i) - beta(4) - beta(5) - - SELECT CASE(TRIM(nonsingularNormal)) - CASE("x") - beta(2) = ( (v(i)-vn*ny)*bsasdq(1,i) + nx*ny*bsasdq(mu,i) & - + (ny**2-1)*bsasdq(mv,i) + ny*nz*bsasdq(mw,i) )/nx - beta(3) = (-(w(i)-vn*nz)*bsasdq(1,i) - nx*nz*bsasdq(mu,i) & - - ny*nz*bsasdq(mv,i) - (nz**2-1)*bsasdq(mw,i) )/nx - CASE("y") - beta(2) = (-(u(i)-vn*nx)*bsasdq(1,i) - (nx**2-1)*bsasdq(mu,i) & - - ny*nx*bsasdq(mv,i) - nx*nz*bsasdq(mw,i) )/ny - beta(3) = ( (w(i)-vn*nz)*bsasdq(1,i) + nx*nz*bsasdq(mu,i) & - + ny*nz*bsasdq(mv,i) + (nz**2-1)*bsasdq(mw,i) )/ny - CASE("z") - beta(2) = ( (u(i)-vn*nx)*bsasdq(1,i) + (nx**2-1)*bsasdq(mu,i) & - + nx*ny*bsasdq(mv,i) + nz*nx*bsasdq(mw,i) )/nz - beta(3) = (-(v(i)-vn*ny)*bsasdq(1,i) - nx*ny*bsasdq(mu,i) & - - (ny**2-1)*bsasdq(mv,i) - nz*ny*bsasdq(mw,i) )/nz - END SELECT - - waveb(1,1) = beta(1) - waveb(mu,1) = beta(1)*(u(i)-a(i)*nx) - waveb(mv,1) = beta(1)*(v(i)-a(i)*ny) - waveb(mw,1) = beta(1)*(w(i)-a(i)*nz) - waveb(5,1) = beta(1)*(enth(i) - vn*a(i)) - sb(1) = vn - a(i) - - waveb(1,2) = beta(4) - SELECT CASE(TRIM(nonsingularNormal)) - CASE("x") - waveb(mu,2) = beta(4)*u(i) + beta(2)*ny - beta(3)*nz - waveb(mv,2) = beta(4)*v(i) - beta(2)*nx - waveb(mw,2) = beta(4)*w(i) + beta(3)*nx - waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & - + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(w(i)*nx-u(i)*nz) - CASE("y") - waveb(mu,2) = beta(4)*u(i) + beta(2)*ny - waveb(mv,2) = beta(4)*v(i) - beta(2)*nx + beta(3)*nz - waveb(mw,2) = beta(4)*w(i) - beta(3)*ny - waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & - + beta(2)*(u(i)*ny-v(i)*nx) + beta(3)*(v(i)*nz-w(i)*ny) - CASE("z") - waveb(mu,2) = beta(4)*u(i) - beta(2)*nz - waveb(mv,2) = beta(4)*v(i) + beta(3)*nz - waveb(mw,2) = beta(4)*w(i) + beta(2)*nx - beta(3)*ny - waveb(5,2) = beta(4)*0.5d0*u2v2w2(i) & - + beta(2)*(w(i)*nx-u(i)*nz) + beta(3)*(v(i)*nz-w(i)*ny) - END SELECT - sb(2) = vn - - waveb(1,3) = beta(5) - waveb(mu,3) = beta(5)*(u(i) + a(i)*nx) - waveb(mv,3) = beta(5)*(v(i) + a(i)*ny) - waveb(mw,3) = beta(5)*(w(i) + a(i)*nz) - waveb(5,3) = beta(5)*(enth(i) + vn*a(i)) - sb(3) = vn + a(i) - - cmbsasdq(:,i) = 0.d0 - cpbsasdq(:,i) = 0.d0 - DO mws = 1,mwaves - IF (sb(mws) < 0.d0)THEN - waveb(:,mws) = waveb(:,mws)*areaRatioMinus - cmbsasdq(:,i) = cmbsasdq(:,i) + sb(mws)*waveb(:,mws) - ELSE - waveb(:,mws) = waveb(:,mws)*areaRatioPlus - cpbsasdq(:,i) = cpbsasdq(:,i) + sb(mws)*waveb(:,mws) - END IF - END DO - END DO - -END SUBROUTINE rptt3 From e1f440e81176a249a50f636521f0f79d9b392f2e Mon Sep 17 00:00:00 2001 From: Carlos Munoz Moncayo Date: Thu, 4 Jul 2024 16:50:44 +0300 Subject: [PATCH 5/7] Update imports --- examples/euler_gravity_3d/rising_hot_sphere.py | 6 +++--- examples/euler_gravity_3d/rising_hot_sphere_spherical.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/euler_gravity_3d/rising_hot_sphere.py b/examples/euler_gravity_3d/rising_hot_sphere.py index b95d006f6..d545e89a2 100755 --- a/examples/euler_gravity_3d/rising_hot_sphere.py +++ b/examples/euler_gravity_3d/rising_hot_sphere.py @@ -11,7 +11,7 @@ density (rho), x,y, and z momentum (rho*u,rho*v,rho*w), and energy. """ import numpy as np -from clawpack.pyclaw.examples.euler_gravity_3d.mappedGrid import euler3d_mappedgrid as mg +from clawpack.riemann.mappedGrid import euler3d_mappedgrid as mg try: from mpi4py import MPI @@ -176,8 +176,8 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ import logging solver.logger.setLevel(logging.DEBUG) - from clawpack.pyclaw.examples.euler_gravity_3d import euler_3d_gmap - solver.rp = euler_3d_gmap + from clawpack import riemann + solver.rp = riemann.euler_3d_gmap solver.num_eqn = 5 solver.num_waves = 3 solver.cfl_max = 0.6 diff --git a/examples/euler_gravity_3d/rising_hot_sphere_spherical.py b/examples/euler_gravity_3d/rising_hot_sphere_spherical.py index 58e62a335..79bd3f9d7 100755 --- a/examples/euler_gravity_3d/rising_hot_sphere_spherical.py +++ b/examples/euler_gravity_3d/rising_hot_sphere_spherical.py @@ -11,7 +11,7 @@ density (rho), x,y, and z momentum (rho*u,rho*v,rho*w), and energy. """ import numpy as np -from clawpack.pyclaw.examples.euler_gravity_3d.mappedGrid import euler3d_mappedgrid as mg +from clawpack.riemann.mappedGrid import euler3d_mappedgrid as mg # Test for MPI, and set sizes accordingly try: @@ -173,8 +173,8 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ import logging solver.logger.setLevel(logging.DEBUG) - from clawpack.pyclaw.examples.euler_gravity_3d import euler_3d_gmap - solver.rp = euler_3d_gmap + from clawpack import riemann + solver.rp = riemann.euler_3d_gmap solver.num_eqn = 5 solver.num_waves = 3 solver.cfl_max = 0.6 From fd2594cd0648d41e8145e6851395ec8881d3d9a5 Mon Sep 17 00:00:00 2001 From: Carlos Munoz Moncayo Date: Thu, 4 Jul 2024 16:55:06 +0300 Subject: [PATCH 6/7] Update meson.build --- src/pyclaw/meson.build | 46 ------------------------------------------ 1 file changed, 46 deletions(-) diff --git a/src/pyclaw/meson.build b/src/pyclaw/meson.build index 3f7f3c038..d7381d885 100644 --- a/src/pyclaw/meson.build +++ b/src/pyclaw/meson.build @@ -326,49 +326,3 @@ py.extension_module( subdir: 'clawpack/pyclaw/examples/advection_reaction_2d', install : true ) - -# 3D Euler with gravity - -ext_name = 'euler_3d_gmap' -ext_srcs = [ - '../../examples/euler_gravity_3d/rpn3_euler_mapgrid.f90', - '../../examples/euler_gravity_3d/rpt3_euler_mapgrid.f90', - '../../examples/euler_gravity_3d/rptt3_euler_mapgrid.f90', -] -f2py_srcs = custom_target( - 'f2py_3d_gmap', - command: [f2py, ext_name], - input: ext_srcs, - output: [ext_name + 'module.c', ext_name + '-f2pywrappers.f'], - build_by_default: true, -) - -py.extension_module( - ext_name, [ext_srcs, f2py_srcs], - incdir_f2py / 'fortranobject.c', - include_directories: inc_np, - dependencies : py_dep, - subdir: 'clawpack/pyclaw/examples/euler_gravity_3d', - install : true -) - -ext_name = 'mappedGrid' -ext_srcs = [ - '../../examples/euler_gravity_3d/euler3d_mappedGrid.f90', -] -f2py_srcs = custom_target( - 'f2py_mappedGrid', - command: [f2py, ext_name], - input: ext_srcs, - output: [ext_name + 'module.c', ext_name + '-f2pywrappers2.f90'], - build_by_default: true, -) - -py.extension_module( - ext_name, [ext_srcs, f2py_srcs], - incdir_f2py / 'fortranobject.c', - include_directories: inc_np, - dependencies : py_dep, - subdir: 'clawpack/pyclaw/examples/euler_gravity_3d', - install : true -) \ No newline at end of file From d84f2ff7bb82f50d50864f3373f2a6fd67d90f28 Mon Sep 17 00:00:00 2001 From: Carlos Munoz Moncayo Date: Fri, 5 Jul 2024 23:41:37 +0300 Subject: [PATCH 7/7] Fix imports and do not keep copy in mem by default --- examples/euler_gravity_3d/rising_hot_sphere.py | 6 +++--- examples/euler_gravity_3d/rising_hot_sphere_spherical.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/euler_gravity_3d/rising_hot_sphere.py b/examples/euler_gravity_3d/rising_hot_sphere.py index d545e89a2..4c1f8942d 100755 --- a/examples/euler_gravity_3d/rising_hot_sphere.py +++ b/examples/euler_gravity_3d/rising_hot_sphere.py @@ -155,7 +155,7 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ use_petsc=False,outdir='./_output',\ output_format='hdf5',file_prefix='equil',disable_output=False,\ mx=mxyz[0],my=mxyz[1],mz=mxyz[2],\ - tfinal=64.0,num_output_times=1): + tfinal=64.0,num_output_times=1,keep_copy=False): if use_petsc: import clawpack.petclaw as pyclaw @@ -177,7 +177,7 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ import logging solver.logger.setLevel(logging.DEBUG) from clawpack import riemann - solver.rp = riemann.euler_3d_gmap + solver.rp = riemann.euler_mapgrid_3D solver.num_eqn = 5 solver.num_waves = 3 solver.cfl_max = 0.6 @@ -340,7 +340,7 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ claw.solver = solver claw.output_format = output_format claw.output_file_prefix = file_prefix - claw.keep_copy = True + claw.keep_copy = keep_copy if disable_output: claw.output_format = None claw.tfinal = tfinal diff --git a/examples/euler_gravity_3d/rising_hot_sphere_spherical.py b/examples/euler_gravity_3d/rising_hot_sphere_spherical.py index 79bd3f9d7..b9604a6e6 100755 --- a/examples/euler_gravity_3d/rising_hot_sphere_spherical.py +++ b/examples/euler_gravity_3d/rising_hot_sphere_spherical.py @@ -174,7 +174,7 @@ def euler3d(kernel_language='Fortran',solver_type='classic',\ import logging solver.logger.setLevel(logging.DEBUG) from clawpack import riemann - solver.rp = riemann.euler_3d_gmap + solver.rp = riemann.euler_mapgrid_3D solver.num_eqn = 5 solver.num_waves = 3 solver.cfl_max = 0.6