big sphere advection, tc1 ok.

Advection2.f90

@@ -383,14 +383,14 @@ subroutine AdvectionRK4(self,aMesh, dt, t, procRank, nProcs, velocityFunction)
 
 ! ! STAGE 1 ONLY : Store velocity and kinetic energy
  do j=1,aParticles%N
- aParticles%ke(j) = sum(self%particlesStage1(:,j)*self%particlesStage1(:,j))
+ !aParticles%ke(j) = sum(self%particlesStage1(:,j)*self%particlesStage1(:,j))
  aParticles%u(:,j) = self%particlesStage1(:,j)
  enddo
  do j=1,self%activePanels%N
- self%activePanels%ke(j) = sum(self%activePanelsStage1(:,j)*self%activePanelsStage1(:,j))
+ !self%activePanels%ke(j) = sum(self%activePanelsStage1(:,j)*self%activePanelsStage1(:,j))
  self%activePanels%u(:,j) = self%activePanelsStage1(:,j)
  enddo
- aMesh%totalKE = 0.5*sum(self%activePanels%ke*self%area)
+ !aMesh%totalKE = 0.5*sum(self%activePanels%ke*self%area)
 
  self%particlesStage1 = dt*self%particlesStage1
  self%passivePanelsStage1 = dt*self%passivePanelsStage1
@@ -664,8 +664,7 @@ function TestCase1Velocity(xyz,t)
 ! solid-body rotation wind field from Williamson, Drake, Hack, Jakob, and Swarztrauber, JCP, 1992
  real(kreal), intent(in) :: xyz(3), t
  real(kreal) :: TestCase1Velocity(3)
- real(kreal), parameter :: u0 = 2.0_kreal*PI*EARTH_RADIUS/(12.0_kreal*ONE_DAY)
- real(kreal) :: lat, lon, raxis, u, v
+ real(kreal), parameter :: u0 = 2.0_kreal*PI/(12.0_kreal*ONE_DAY)
 
  TestCase1Velocity(1) = -u0*xyz(2)*cos(alpha)
  TestCase1Velocity(2) = u0*(xyz(1)*cos(alpha) - xyz(3)*sin(alpha))

BVEVorticity.f90

@@ -31,6 +31,7 @@ module BVESetupModule
 public New, Delete
 public InitSolidBodyRotation, SetSolidBodyRotationOnMesh
 public SOLID_BODY_NINT, SOLID_BODY_NREAL
+public NullVorticity
 
 !
 !----------------

CubedSphereTest2.f90

@@ -20,6 +20,7 @@ program CubedSphereTest
 use EdgesModule
 use PanelsModule
 use VTKOutputModule
+use AdvectionModule
 
 implicit none
 
@@ -35,6 +36,8 @@ program CubedSphereTest
 type(Logger) :: exeLog
 character(len=56) :: logString
 
+integer(kint) :: j
+
 !
 ! set mesh parameters
 !
@@ -56,6 +59,18 @@ program CubedSphereTest
  call LogStats(sphere,exeLog,trim(logString))
  call LogMessage(exeLog,TRACE_LOGGING_LEVEL,'Surf. area should be = ',4.0_kreal*PI*EARTH_RADIUS*EARTH_RADIUS)
  write(6,'(A,F24.15)') "surf area error = ",abs(4.0_kreal*PI*EARTH_RADIUS*EARTH_RADIUS - sum(spherePanels%area))
+!
+ do j=1,sphereParticles%N
+ sphereParticles%u(:,j) = TestCase1Velocity(sphereParticles%x(:,j),0.0_kreal)
+ enddo
+ do j=1,spherePanels%N
+ if ( .NOT. spherePanels%hasChildren(j) ) then
+ spherePanels%u(:,j) = TestCase1Velocity(spherePanels%x(:,j),0.0_kreal)
+ endif
+ enddo
+!
+
+
 
 !
 ! Define VTKOutput

Makefile

@@ -16,21 +16,30 @@ SHELL = /bin/bash
 #MKL_COMPILE=-openmp -I$(MKLROOT)/include/intel64/lp64 -I$(MKLROOT)/include
 #----------------#
 
+#----------------#
+# vortex.math.lsa.umich.edu
+FF = ifort
+FF_FLAGS = -g -O0 -check bounds -check pointer -check uninit -traceback -warn all -debug extended -openmp
+MKL_ROOT=/usr/local/intel/Compiler/11.1/056/mkl
+MKL_LINK=-L$(MKL_ROOT)/lib $(MKL_ROOT)/lib/libmkl_lapack95_lp64.a -lmkl_intel_lp64 -lmk_intel_thread -lmkl_core -lpthread -lm
+MKL_COMPILE=-openmp -I$(MKL_ROOT)/include/intel64/lp64 -I$(MKL_ROOT)/include
+
+
 #--------------#
 # TANK DESKTOP #
 
-FF=ifort
-FF_FLAGS=-g -traceback -warn all -debug extended
+#FF=ifort
+#FF_FLAGS=-g -traceback -warn all -debug extended
 #FF_FLAGS=-O2 -warn all -opt-report 1
-VTK_INCLUDE=/usr/local/include/vtk-5.10
-VTK_LIB_DIR=/usr/local/lib/vtk-5.10
-VTK_LIBS=-lvtkCommon -lvtkGraphics -lvtkRendering -lvtkViews -lvtkWidgets -lvtkImaging -lvtkHybrid -lvtkIO -lvtkFiltering
-MKLROOT=/opt/intel/mkl
-MKL_THREADING_LAYER=intel
-MKL_INTERFACE_LAYER=lp64
-MKL_LINK=-L$(MKLROOT)/lib $(MKLROOT)/lib/libmkl_lapack95_lp64.a \
--lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread -lm
-MKL_COMPILE=-openmp -I/opt/intel/mkl/include/intel64/lp64 -I/opt/intel/mkl/include 
+#VTK_INCLUDE=/usr/local/include/vtk-5.10
+#VTK_LIB_DIR=/usr/local/lib/vtk-5.10
+#VTK_LIBS=-lvtkCommon -lvtkGraphics -lvtkRendering -lvtkViews -lvtkWidgets -lvtkImaging -lvtkHybrid -lvtkIO -lvtkFiltering
+#MKLROOT=/opt/intel/mkl
+#MKL_THREADING_LAYER=intel
+#MKL_INTERFACE_LAYER=lp64
+#MKL_LINK=-L$(MKLROOT)/lib $(MKLROOT)/lib/libmkl_lapack95_lp64.a \
+#-lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread -lm
+#MKL_COMPILE=-openmp -I/opt/intel/mkl/include/intel64/lp64 -I/opt/intel/mkl/include 
 
 #--------------#
 
@@ -43,10 +52,8 @@ MKL_COMPILE=-openmp -I/opt/intel/mkl/include/intel64/lp64 -I/opt/intel/mkl/inclu
  $(FF) -c $(FF_FLAGS) $<
 %.exe : 
  $(FF) $(FF_FLAGS) -o $@ $< `mpif90 -showme:link` $(MKL_COMPILE) $^
-
 clean:
  rm *.o *.mod *genmod.f90
-
 cleanx:
  rm *.exe 
 
@@ -62,7 +69,7 @@ INTERP_OBJS = $(BASE_OBJS) $(MESH_OBJS) ssrfpack.o stripack.o STRIPACKInterface2
 interp_objs: $(BASE_OBJS) $(MESH_OBJS) ssrfpack.o stripack.o STRIPACKInterface2.o SSRFPACKInterface2.o
 OUTPUT_OBJS = VTKOutput.o $(BASE_OBJS) $(MESH_OBJS)
 TEST_CASE_OBJS = Tracers.o BVEVorticity.o
-ADVECTION_OBJS = $(BASE_OBJS) $(MESH_OBJS) $(INTERP_OBJS) $(OUTPUT_OBJS) $(TEST_CASE_OBJS) Advection2.o
+ADVECTION_OBJS = $(BASE_OBJS) $(MESH_OBJS) $(INTERP_OBJS) $(OUTPUT_OBJS) $(TEST_CASE_OBJS) Advection2.o RefineRemesh2.o
 
 
 #############################################################
@@ -79,7 +86,7 @@ TestCase1.o: TestCase1.f90 $(ADVECTION_OBJS)
 #############################################################
 ## UNIT TEST EXECUTABLES
 #############################################################
-cubedSphereTestSerial.exe: CubedSphereTest2.o $(BASE_OBJS) $(MESH_OBJS) $(INTERP_OBJS) $(OUTPUT_OBJS)
+cubedSphereTestSerial.exe: CubedSphereTest2.o $(BASE_OBJS) $(MESH_OBJS) $(INTERP_OBJS) $(OUTPUT_OBJS) Advection2.o
  $(FF) $(FF_FLAGS) -o $@ $^ `mpif90 -showme:link` $(MKL_COMPILE) 
 
 #############################################################

RefineRemesh2.f90

@@ -107,8 +107,8 @@ subroutine SetVorticityOnMesh(genMesh,genVort)
 !
 subroutine NewPrivate(self, limit, maxTol, varTol, type, tracerID)
  type(RefinementSetup), intent(out) :: self
- real(kreal), intent(in) :: maxTol, varTol, limit
- integer(kint), intent(in) :: type
+ real(kreal), intent(in) :: maxTol, varTol
+ integer(kint), intent(in) :: type, limit
  integer(kint), intent(in), optional :: tracerID
 
  if ( .NOT. logInit ) call InitLogger(log, procRank)
@@ -764,6 +764,11 @@ subroutine FlagPanelsForFlowMapRefinement(refineFlag,aMesh,refineFlowMap,startIn
  enddo
 end subroutine
 
+subroutine FlagPanelsAtRefinementBoundaries(refineFlag,aMesh)
+ logical(klog), intent(inout) :: refineFlag(:)
+ type(SphereMesh), intent(in) :: aMesh
+end subroutine
+
 
 subroutine InitLogger(aLog,rank)
  type(Logger), intent(inout) :: aLog

TestCase1.f90

@@ -24,6 +24,8 @@ program TestCase1
 
 implicit none
 
+include 'mpif.h'
+
 !
 ! mesh variables
 !
@@ -34,7 +36,7 @@ program TestCase1
 !
 type(RK4Data) :: advRK4
 real(kreal) :: t, dt, tfinal
-integer(kint) :: timesteps, timeJ
+integer(kint) :: timesteps, timeJ, timestepsPerPeriod
 !
 ! tracer variables
 !
@@ -54,13 +56,13 @@ program TestCase1
 !
 ! User input
 !
-character(len = 128) :: namelistFile = 'testcase1.namelist'
+character(len = 128) :: namelistInputFile = 'testcase1.namelist'
 integer(kint) :: readStat
 !
 ! Output variables
 !
 type(VTKSource) :: vtkOut
-character(len = 128) :: vtkRoot, vtkFile
+character(len = 128) :: vtkRoot, vtkFile, outputDir, jobPrefix
 character(len = 56 ) :: amrString
 integer(kint) :: frameCounter, frameOut
 !
@@ -74,13 +76,13 @@ program TestCase1
 !
 integer(kint) :: i, j, k
 integer(kint) :: mpiErrCode
-integer(kint), parameter :: BROADCAST_INT_SIZE = 5, BROADCAST_REAL_SIZE = 4
+integer(kint), parameter :: BROADCAST_INT_SIZE = 6, BROADCAST_REAL_SIZE = 3
 integer(kint) :: broadcastIntegers(BROADCAST_INT_SIZE)
 real(kreal) :: broadcastReals(BROADCAST_REAL_SIZE)
 type(BVESetup) :: nullVort
 
 namelist /sphereDefine/ panelKind, initNest, AMR, tracerMaxTol, tracerVarTol, refineMentLimit
-namelist /timeStepping/ tfinal, dt, remeshInterval
+namelist /timeStepping/ tfinal, timestepsPerPeriod,  remeshInterval
 namelist /fileIO/ outputDir, jobPrefix, frameOut
 
 
@@ -102,11 +104,12 @@ program TestCase1
 initNest = broadcastIntegers(3)
 refinementLimit = broadcastIntegers(4)
 remeshInterval = broadcastIntegers(5)
+timestepsPerPeriod = broadCastIntegers(6)
 call MPI_BCAST(broadcastReals,BROADCAST_REAL_SIZE,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpiErrCode)
 tracermaxTol = broadcastReals(1)
 tracerVarTol = broadcastReals(2)
 tfinal = broadcastReals(3)*ONE_DAY ! convert tfinal to seconds
-dt = broadcastreals(4)
+
 !
 ! define test case in accordance with Williamson et al., JCP 1992
 !
@@ -125,8 +128,8 @@ program TestCase1
 call New(sphere,panelKind,initNest,AMR,nTracer,problemKind)
 call SetCosineBellTracerOnMesh(sphere,cosBell)
 if ( AMR > 0 ) then
- call New(tracerRefine,refinementLimit,tracermaxTol,tracerVarTol,REFINE_TRACER,tracerID)
- call InitialRefinement(sphere,tracerRefine,SetCosineBellTracerOnMesh, cosBel, &
+ call New(tracerRefine,refinementLimit,tracermaxTol,tracerVarTol,TRACER_REFINE,tracerID)
+ call InitialRefinement(sphere,tracerRefine,SetCosineBellTracerOnMesh, cosBell, &
  tracerRefine, NullVorticity,nullVort)
  if ( panelKind == QUAD_PANEL) then 
  write(amrString,'(A,I1,A,I0.2)') 'quadAMR_',initNest,'to',initNest+refinementLimit
@@ -154,7 +157,8 @@ program TestCase1
 ! intialize timestepping
 !
 call New(advRK4,sphere,numProcs)
-timesteps = floor(tfinal,dt)
+dt = 12.0_kreal*ONE_DAY/real(timestepsPerPeriod,kreal)
+timesteps = floor(tfinal/dt)
 t = 0.0_kreal
 remeshCounter = 0
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -175,7 +179,7 @@ program TestCase1
  ! build new mesh
  !
  call LagrangianRemesh(sphere, NullVorticity, nullVort, tracerRefine, &
- SetCosineBellOnMesh, cosBell, tracerRefine, &
+ SetCosineBellTracerOnMesh, cosBell, tracerRefine, &
  tracerRefine)
  !
  ! create new associated objects
@@ -189,14 +193,14 @@ program TestCase1
  !
  ! advance timestep
  !
- call AdvectionRK4(advRK4,sphere, dt, t, procRank, numProcs, LauritzenEtAlNonDivergentWind)
+ call AdvectionRK4(advRK4,sphere, dt, t, procRank, numProcs, TestCase1Velocity)
 
  t = real(timeJ+1,kreal)*dt
 
  !
  ! output timestep data
  !
- if ( procRank == 0 .AND. mod(timeJ+1,frameOUt) == 0 ) then
+ if ( procRank == 0 .AND. mod(timeJ+1,frameOUt) == 0  ) then
  call LogMessage(exeLog,TRACE_LOGGING_LEVEL,'day = ',t/ONE_DAY)
  write(vtkFile,'(A,I0.4,A)') trim(vtkRoot), frameCounter, '.vtk'
  call UpdateFileName(vtkOut,vtkFile)
@@ -219,7 +223,7 @@ program TestCase1
 call Delete(sphere)
 call Delete(cosBell)
 call Delete(exeLog)
-
+call MPI_FINALIZE(mpiErrCode)
 contains
 
 subroutine ReadNamelistfile(rank)
@@ -240,11 +244,11 @@ subroutine ReadNamelistfile(rank)
  broadcastIntegers(3) = initNest
  broadcastIntegers(4) = refinementLimit
  broadcastIntegers(5) = remeshInterval
-
+ broadcastIntegers(6) = timestepsPerPeriod
+
  broadcastReals(1) = tracerMaxTol
  broadcastReals(2) = tracerVarTol
- broadcastReals(3) = tfinal
- broadcastReals(4) = dt 
+ broadcastReals(3) = tfinal 
  endif
 end subroutine
 
@@ -256,7 +260,7 @@ subroutine InitLogger(alog,rank)
  else
  call New(aLog,WARNING_LOGGING_LEVEL)
  endif
-
+ write(logKey,'(A,I0.2,A)') 'EXE_LOG',procRank,' : '
 end subroutine
 
-end program
+end program

Tracers.f90

@@ -77,9 +77,9 @@ module TracerSetupModule
 ! Standard methods : Constructor / Destructor
 !----------------
 !
-subroutine NewPrivate(self, nInt, nReal, tracerID)
+subroutine NewPrivate(self, nInt, nReal)
  type(TracerSetup), intent(out) :: self
- integer(kint), intent(in) :: nINt, nReal, tracerID
+ integer(kint), intent(in) :: nINt, nReal
 
  if (.NOT. logInit) call InitLogger(log,procRank)
 
@@ -94,8 +94,7 @@ subroutine NewPrivate(self, nInt, nReal, tracerID)
  self%reals = 0.0_kreal
  else
  nullify(self%reals)
- endif 
- self%tracerID = tracerID 
+ endif 
 end subroutine
 
 subroutine DeletePrivate(self)

VTKOutput.f90

@@ -171,11 +171,11 @@ subroutine vtkOutput(self,aMesh)
  !
  write(WRITE_UNIT_1,'(A,I8,A)') 'POINTS ',self%nPoints,' double'
  do j=1,aparticles%N
- write(WRITE_UNIT_1,'(3F24.15)') aParticles%x(1,j),aParticles%x(2,j),aParticles%x(3,j)
+ write(WRITE_UNIT_1,'(3F24.8)') aParticles%x(1,j),aParticles%x(2,j),aParticles%x(3,j)
  enddo
  do j=1,aPanels%N
  if ( .NOT. aPanels%hasChildren(j) ) then
- write(WRITE_UNIT_1,'(3F24.15)') aPanels%x(1,j), aPanels%x(2,j), aPanels%x(3,j)
+ write(WRITE_UNIT_1,'(3F24.8)') aPanels%x(1,j), aPanels%x(2,j), aPanels%x(3,j)
  endif
  enddo
 
@@ -202,11 +202,11 @@ subroutine vtkOutput(self,aMesh)
  write(WRITE_UNIT_1,'(A)') trim(dataString)
  write(WRITE_UNIT_1,'(A)') 'LOOKUP_TABLE default'
  do j=1,aParticles%N
- write(WRITE_UNIT_1,'(3F24.15)') aParticles%x0(:,j)
+ write(WRITE_UNIT_1,'(3F24.8)') aParticles%x0(:,j)
  enddo
  do j=1,aPanels%N
  if (.NOT. aPanels%hasChildren(j) ) then
- write(WRITE_UNIT_1,'(3F24.15)') aPanels%x0(:,j)
+ write(WRITE_UNIT_1,'(3F24.8)') aPanels%x0(:,j)
  endif
  enddo
 
@@ -263,7 +263,7 @@ subroutine vtkOutput(self,aMesh)
  enddo
  do j=1,aPanels%N
  if ( .NOT. aPanels%hasChildren(j) ) then
- write(WRITE_UNIT_1,'(3F24.15)') aPanels%u(:,j)
+ write(WRITE_UNIT_1,'(3F24.8)') aPanels%u(:,j)
  endif
  enddo
  endif

testcase1.namelist

@@ -9,12 +9,12 @@
 
 & timeStepping
  tfinal = 12 ! days
- dt = 2592   ! seconds
- remeshInterval = 20
+ timestepsPerPeriod = 400 
+ remeshInterval = 30
 /
 
 & fileIO
- outputDir = '~/Desktop/testCase1'
- jobPrefix = 'tc1_adv_rm20'
+ outputDir = '~/BigSphere/tc1'
+ jobPrefix = 'tc1_adv_rm30'
  frameOut = 4
 /