Skip to content

Commit

Permalink
gaussvort, almost. not sure about units.
Browse files Browse the repository at this point in the history
  • Loading branch information
pbosler committed Jan 27, 2014
1 parent cd7f23f commit 8a42e24
Show file tree
Hide file tree
Showing 7 changed files with 105 additions and 91 deletions.
8 changes: 8 additions & 0 deletions BVEDirectSum.f90
Original file line number Diff line number Diff line change
Expand Up @@ -30,12 +30,15 @@ module BVEDirectSumModule
public BVERK4Data
public New, Delete
public BVERK4Timestep
public VELOCITY_SMOOTH

!
!----------------
! Types and module constants
!----------------
!
real(kreal), protected, save :: VELOCITY_SMOOTH = 0.01_kreal

type BVERK4Data
! MPI load balancing
integer(kint), pointer :: particlesIndexStart(:) => null(), &
Expand Down Expand Up @@ -294,6 +297,11 @@ subroutine DeletePrivate(self)
! Public member functions
!----------------
!
subroutine SetVelocitySmoothingParameter(newSmooth)
real(kreal), intent(in) :: newSmooth
VELOCITY_SMOOTH = newSmooth
end subroutine

subroutine BVERK4Timestep(self,aMesh,dt,procRank, nProcs)
! Advances aMesh one time step forward using 4th order Runge-Kutta
!
Expand Down
12 changes: 6 additions & 6 deletions BVESingleGaussianVortex.f90
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ program BVESolidBody

namelist /sphereDefine/ panelKind, initNest, AMR, tracerMaxTol, tracerVarTol, &
refineMentLimit, circMaxTol, vortVarTol, lagVarTol
namelist /vorticityDefine/ vortLat, vortLon, bb, maxVort
namelist /vorticityDefine/ vortLat, vortLon, bb, maxVort
namelist /timeStepping/ tfinal, dt, remeshInterval
namelist /fileIO/ outputDir, jobPrefix, frameOut

Expand Down Expand Up @@ -127,7 +127,7 @@ program BVESolidBody
vortLat = broadcastReals(8)
vortLon = broadcastReals(9)
bb = broadcastReals(10)
maxVort = broadcastReals(11)
maxVort = broadcastReals(11)*OMEGA


!
Expand Down Expand Up @@ -159,7 +159,7 @@ program BVESolidBody
call New(flowMapREfine,refinementLimit,100000.0_kreal,lagVarTol,FLOWMAP_REFINE)
call InitialRefinement(sphere,tracerRefine,SetCosineBellTracerOnMesh, cosBell, &
vortRefine, SetSingleGaussianVortexOnMesh,gaussVort)
call SetInitialLatitudeTracerOnMesh(sphere,2)
call SetInitialLatitudeTracerOnMesh(sphere,2)
if ( panelKind == QUAD_PANEL) then
write(amrString,'(A,I1,A,I0.2)') 'quadAMR_',initNest,'to',initNest+refinementLimit
endif
Expand Down Expand Up @@ -241,7 +241,7 @@ program BVESolidBody
! output final data
!

! TO DO : output final data
! TO DO : output final data



Expand Down Expand Up @@ -269,8 +269,8 @@ subroutine ReadNamelistfile(rank)
endif
read(READ_UNIT,nml=sphereDefine)
rewind(READ_UNIT)
read(READ_UNIT,nml=vorticityDefine)
rewind(READ_UNIT)
read(READ_UNIT,nml=vorticityDefine)
rewind(READ_UNIT)
read(READ_UNIT,nml=timestepping)
rewind(READ_UNIT)
read(READ_UNIT,nml=fileIO)
Expand Down
26 changes: 11 additions & 15 deletions BVEVorticity.f90
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ module BVESetupModule
public NullVorticity
public InitSingleGaussianVortex, SetSingleGaussianVortexOnMesh
public SINGLE_VORTEX_NINT, SINGLE_VORTEX_NREAL
public GAUSS_CONST

!
!----------------
Expand All @@ -45,6 +46,7 @@ module BVESetupModule
real(kreal), pointer :: reals(:) => null() ! real numbers required by test case definitions
end type

real(kreal), protected, save :: GAUSS_CONST = 0.0_kreal
integer(kint), parameter :: SOLID_BODY_NINT = 0, SOLID_BODY_NREAL = 1
integer(kint), parameter :: SINGLE_VORTEX_NINT = 0, SINGLE_VORTEX_NREAL = 4

Expand Down Expand Up @@ -173,14 +175,13 @@ subroutine SetSingleGaussianVortexOnMesh(aMesh,gaussVort)
type(Particles), pointer :: aParticles
type(Panels), pointer :: aPanels
real(kreal), allocatable :: gVort(:)

aParticles => aMesh%particles
aPanels => aMesh%panels

xyzCent = EARTH_RADIUS*[ cos(gaussVort%reals(1))*cos(gaussVort%reals(2)), &
cos(gaussVort%reals(1))*sin(gaussVort%reals(2)), &
sin(gaussVort%reals(1)) ]
print *, xyzCent
!
! set Gaussian constant if not set already
!
Expand All @@ -193,22 +194,17 @@ subroutine SetSingleGaussianVortexOnMesh(aMesh,gaussVort)
gaussVort%reals(3), gaussVort%reals(4))
endif
enddo
GAUSS_CONST = sum(gVort*aPanels%area(1:aPanels%N))/(4.0_kreal*PI*EARTH_RADIUS)

call LogMessage(log,DEBUG_LOGGING_LEVEL,'GAUSS_CONST = ',GAUSS_CONST)
call LogMessage(log,DEBUG_LOGGING_LEVEL,'n low level panels = ',aPanels%N - count(aPanels%hasChildren))
call LogMessage(log,DEBUG_LOGGING_LEVEL,'gauss int = ',sum(gVort*aPanels%area(1:aPanels%N)))

GAUSS_CONST = sum(gVort*aPanels%area(1:aPanels%N))/(4.0_kreal*PI*EARTH_RADIUS*EARTH_RADIUS)
deallocate(gVort)
endif

do j=1,aParticles%N
aParticles%absVort(j) = GeneralGaussian(aParticles%x0(:,j), xyzCent, &
gaussVort%reals(3), gaussVort%reals(4)) &
- GAUSS_CONST + 2.0_kreal*OMEGA*aParticles%x0(3,j)/EARTH_RADIUS
aParticles%relVort(j) = aParticles%absVort(j) - 2.0_kreal*OMEGA*aParticles%x(3,j)/EARTH_RADIUS
aParticles%relVort(j) = aParticles%absVort(j) - 2.0_kreal*OMEGA*aParticles%x(3,j)/EARTH_RADIUS
enddo

do j=1,aPanels%N
if ( aPanels%hasChildren(j) ) then
aPanels%absVort(j) = 0.0_kreal
Expand All @@ -217,9 +213,9 @@ subroutine SetSingleGaussianVortexOnMesh(aMesh,gaussVort)
aPanels%absVort(j) = GeneralGaussian(aPanels%x0(:,j), xyzCent, &
gaussVort%reals(3), gaussVort%reals(4)) &
- GAUSS_CONST + 2.0_kreal*OMEGA*aPanels%x0(3,j)/EARTH_RADIUS
aPanels%relVort(j) = aPanels%absVort(j) - 2.0_kreal*OMEGA*aPanels%x(3,j)/EARTH_RADIUS
aPanels%relVort(j) = aPanels%absVort(j) - 2.0_kreal*OMEGA*aPanels%x(3,j)/EARTH_RADIUS
endif
enddo
enddo
end subroutine

subroutine NullVorticity(aMesh,nullVort)
Expand All @@ -244,7 +240,7 @@ function SolidBodyX(xyz,rotationRate)
function GeneralGaussian(xyz,xyzCent, bb, maxVal)
real(kreal) :: GeneralGaussian
real(kreal), intent(in) :: xyz(3), xyzCent(3), bb, maxVal
GeneralGaussian = maxVal * exp( -2.0_kreal*bb*bb*( EARTH_RADIUS*EARTH_RADIUS - sum(xyz*xyzCent)))
GeneralGaussian = maxVal * exp( -2.0_kreal*bb*bb/EARTH_RADIUS/EARTH_RADIUS *( sum((xyz-xyzCent)*(xyz-xyzCent))))
end function

subroutine InitLogger(aLog,rank)
Expand Down
3 changes: 0 additions & 3 deletions NumberKinds3.f90
Original file line number Diff line number Diff line change
Expand Up @@ -49,9 +49,6 @@ module NumberKindsModule
integer(KINT), save :: numProcs = 1, &
procRank = 0

! Test case variables
real(kreal), save :: GAUSS_CONST = 0.0_kreal,&
VELOCITY_SMOOTH = 0.01_kreal

end module NumberKindsModule

Loading

0 comments on commit 8a42e24

Please sign in to comment.