Spherical

Makefile

@@ -339,8 +339,11 @@ Write_netcdf.o:./src/Write_netcdf.cu
 
 Write_txtlog.o:./src/Write_txtlog.cpp
 	$(EXEC) $(NVCC) $(INCLUDES) $(ALL_CCFLAGS) $(GENCODE_FLAGS) -o $@ -c $<
+
+Spherical.o:./src/Spherical.cu
+	$(EXEC) $(NVCC) $(INCLUDES) $(ALL_CCFLAGS) $(GENCODE_FLAGS) -o $@ -c $<
 
-BG_Flood: BG_Flood.o AdaptCriteria.o Adaptation.o Advection.o Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o Write_txtlog.o Poly.o
+BG_Flood: BG_Flood.o AdaptCriteria.o Adaptation.o Advection.o Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o Write_txtlog.o Poly.o Spherical.o
 	$(EXEC) $(NVCC) $(ALL_LDFLAGS) $(GENCODE_FLAGS) -o $@ $+ $(LIBRARIES)
 	$(EXEC) mkdir -p ../../bin/$(TARGET_ARCH)/$(TARGET_OS)/$(BUILD_TYPE)
 	$(EXEC) cp $@ ../../bin/$(TARGET_ARCH)/$(TARGET_OS)/$(BUILD_TYPE)
@@ -349,7 +352,7 @@ run: build
 	$(EXEC) ./BG_Flood
 
 clean:
-	rm -f BG_Flood AdaptCriteria.o Adaptation.o Advection.o BG_Flood.o Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o Write_txtlog.o Poly.o
+	rm -f BG_Flood AdaptCriteria.o Adaptation.o Advection.o BG_Flood.o Boundary.o ConserveElevation.o FlowCPU.o FlowGPU.o Friction.o Gradients.o GridManip.o Halo.o InitEvolv.o InitialConditions.o Kurganov.o Mainloop.o Meanmax.o MemManagement.o Mesh.o ReadForcing.o ReadInput.o Read_netcdf.o Reimann.o Setup_GPU.o Testing.o Updateforcing.o Util_CPU.o Write_netcdf.o Write_txtlog.o Poly.o Spherical.o
 	rm -rf ../../bin/$(TARGET_ARCH)/$(TARGET_OS)/$(BUILD_TYPE)/template
 
 clobber: clean

doc/codemap.md

@@ -0,0 +1,7 @@
+```mermaid
+  graph TD;
+      A-->B;
+      A-->C;
+      B-->D;
+      C-->D;
+```

src/Adaptation.h

@@ -13,6 +13,7 @@
 #include "InitialConditions.h"
 #include "Testing.h"
 
+
 template <class T> void Adaptation(Param& XParam, Forcing<float> XForcing, Model<T>& XModel);
 template <class T> void InitialAdaptation(Param& XParam, Forcing<float> &XForcing, Model<T>& XModel);
 template <class T> bool refinesanitycheck(Param XParam, BlockP<T> XBlock, bool*& refine, bool*& coarsen);

src/Advection.cu

@@ -37,34 +37,58 @@ struct SharedMemory<double>
 
 template <class T>__global__ void updateEVGPU(Param XParam, BlockP<T> XBlock, EvolvingP<T> XEv, FluxP<T> XFlux, AdvanceP<T> XAdv)
 {
-	
-	unsigned int halowidth = XParam.halowidth;
-	unsigned int blkmemwidth = blockDim.x + halowidth * 2;
+
+	int halowidth = XParam.halowidth;
+	int blkmemwidth = blockDim.x + halowidth * 2;
 	//unsigned int blksize = blkmemwidth * blkmemwidth;
-	unsigned int ix = threadIdx.x;
-	unsigned int iy = threadIdx.y;
-	unsigned int ibl = blockIdx.x;
-	unsigned int ib = XBlock.active[ibl];
+	int ix = threadIdx.x;
+	int iy = threadIdx.y;
+	int ibl = blockIdx.x;
+	int ib = XBlock.active[ibl];
+
+	int lev = XBlock.level[ib];
 
 	//T eps = T(XParam.eps);
-	T delta = calcres(T(XParam.dx), XBlock.level[ib]);
+	T delta = calcres(T(XParam.delta), lev);
 	T g = T(XParam.g);
-
 
-	T fc = (T)XParam.lat * pi / T(21600.0);
+	T ybo = T(XParam.yo + XBlock.yo[ib]);
+
+
+	T fc = 0.0;// XParam.spherical ? sin((ybo + calcres(T(XParam.dx), lev) * iy) * pi / 180.0) * pi / T(21600.0) : sin(T(XParam.lat * pi / 180.0)) * pi / T(21600.0); // 2*(2*pi/24/3600)
+	// fc should be pi / T(21600.0) * sin(phi)
+
 
 	int iright, itop;
-	
+
 	int i = memloc(halowidth, blkmemwidth, ix, iy, ib);
-	
+
 	iright = memloc(halowidth, blkmemwidth, ix + 1, iy, ib);
 	itop = memloc(halowidth, blkmemwidth, ix, iy + 1, ib);
 
-
+	T yup = T(iy) + T(1.0);
+
+	if (iy == XParam.blkwidth - 1)
+	{
+		if (XBlock.level[XBlock.TopLeft[ib]] > XBlock.level[ib])
+		{
+			yup = iy + 0.75;
+		}
+		if (XBlock.level[XBlock.TopLeft[ib]] < XBlock.level[ib])
+		{
+			yup = iy + 1.000;
+		}
+	}
+
 
-	T cm = T(1.0);// 0.1;
+
+	T cm = XParam.spherical ? calcCM(T(XParam.Radius), delta, ybo, iy) : T(1.0);
 	T fmu = T(1.0);
-	T fmv = T(1.0);
+	T fmv = XParam.spherical ? calcFM(T(XParam.Radius), delta, ybo, T(iy)) : T(1.0);
+	T fmup = T(1.0);
+	T fmvp = XParam.spherical ? calcFM(T(XParam.Radius), delta, ybo, yup) : T(1.0);
+
+
 
 	T hi = XEv.h[i];
 	T uui = XEv.u[i];
@@ -73,24 +97,33 @@ template <class T>__global__ void updateEVGPU(Param XParam, BlockP<T> XBlock, Ev
 
 	T cmdinv, ga;
 
-	cmdinv = T(1.0)/ (cm * delta);
+	cmdinv = T(1.0) / (cm * delta);
 	ga = T(0.5) * g;
-		
+
 
 	XAdv.dh[i] = T(-1.0) * (XFlux.Fhu[iright] - XFlux.Fhu[i] + XFlux.Fhv[itop] - XFlux.Fhv[i]) * cmdinv;
-		
+
 
 
 	//double dmdl = (fmu[xplus + iy*nx] - fmu[i]) / (cm * delta);
 	//double dmdt = (fmv[ix + yplus*nx] - fmv[i]) / (cm  * delta);
-	T dmdl = (fmu - fmu) / (cm * delta);// absurd if not spherical!
-	T dmdt = (fmv - fmv) / (cm * delta);
+	T dmdl = (fmup - fmu) * cmdinv;// absurd if not spherical!
+	T dmdt = (fmvp - fmv) * cmdinv;;
 	T fG = vvi * dmdl - uui * dmdt;
 	XAdv.dhu[i] = (XFlux.Fqux[i] + XFlux.Fquy[i] - XFlux.Su[iright] - XFlux.Fquy[itop]) * cmdinv + fc * hi * vvi;
 	XAdv.dhv[i] = (XFlux.Fqvy[i] + XFlux.Fqvx[i] - XFlux.Sv[itop] - XFlux.Fqvx[iright]) * cmdinv - fc * hi * uui;
-		
+
 	XAdv.dhu[i] += hi * (ga * hi * dmdl + fG * vvi);// This term is == 0 so should be commented here
 	XAdv.dhv[i] += hi * (ga * hi * dmdt - fG * uui);// Need double checking before doing that
+
+	if (XBlock.level[XBlock.TopLeft[ib]] < XBlock.level[ib])
+	{
+		if (iy == XParam.blkwidth - 1)
+		{
+			printf("XAdv.dh[i]=%e, XAdv.dhv[i]=%e, XAdv.dhv[i]=%e, dmdt=%e \n", XAdv.dh[i], XAdv.dhv[i], XAdv.dhu[i], dmdt);
+		}
+
+	}
 
 }
 template __global__ void updateEVGPU<float>(Param XParam, BlockP<float> XBlock, EvolvingP<float> XEv, FluxP<float> XFlux, AdvanceP<float> XAdv);
@@ -103,22 +136,28 @@ template <class T>__host__ void updateEVCPU(Param XParam, BlockP<T> XBlock, Evol
 	//T eps = T(XParam.eps);
 	T delta;
 	T g = T(XParam.g);
+
+	T ybo;
 
 
-	int ib;
+	int ib,lev;
 	int halowidth = XParam.halowidth;
 	int blkmemwidth = XParam.blkmemwidth;
 
 	for (int ibl = 0; ibl < XParam.nblk; ibl++)
 	{
 		ib = XBlock.active[ibl];
-		delta = calcres(T(XParam.dx), XBlock.level[ib]);
+		lev = XBlock.level[ib];
+		delta = calcres(T(XParam.delta), lev);
+
+		ybo = XParam.yo + XBlock.yo[ib];
+
 		for (int iy = 0; iy < XParam.blkwidth; iy++)
 		{
 			for (int ix = 0; ix < XParam.blkwidth; ix++)
 			{
 
-				T fc = T(XParam.lat * pi / 21600.0);
+				T fc = XParam.spherical ? sin((ybo + calcres(T(XParam.dx), lev) * iy) * pi / 180.0) * pi / T(21600.0) : sin(T(XParam.lat * pi / 180.0)) * pi / T(21600.0); // 2*(2*pi/24/3600)
 
 				int iright, itop;
 
@@ -128,10 +167,27 @@ template <class T>__host__ void updateEVCPU(Param XParam, BlockP<T> XBlock, Evol
 				itop = memloc(halowidth, blkmemwidth, ix, iy + 1, ib);
 
 
+				T yup = T(iy) + T(1.0);
+
+				if (iy == XParam.blkwidth - 1)
+				{
+					if (XBlock.level[XBlock.TopLeft[ib]] > XBlock.level[ib])
+					{
+						yup = iy + 0.75;
+					}
+					if (XBlock.level[XBlock.TopLeft[ib]] < XBlock.level[ib])
+					{
+						yup = iy + 1.5;
+					}
+				}
+
+
 
-				T cm = T(1.0);// 0.1;
+				T cm = XParam.spherical ? calcCM(T(XParam.Radius), delta, ybo, iy) : T(1.0);
 				T fmu = T(1.0);
-				T fmv = T(1.0);
+				T fmv = XParam.spherical ? calcFM(T(XParam.Radius), delta, ybo, T(iy)) : T(1.0);
+				T fmup = T(1.0);
+				T fmvp = XParam.spherical ? calcFM(T(XParam.Radius), delta, ybo, yup) : T(1.0);
 
 				T hi = XEv.h[i];
 				T uui = XEv.u[i];
@@ -150,8 +206,8 @@ template <class T>__host__ void updateEVCPU(Param XParam, BlockP<T> XBlock, Evol
 
 				//double dmdl = (fmu[xplus + iy*nx] - fmu[i]) / (cm * delta);
 				//double dmdt = (fmv[ix + yplus*nx] - fmv[i]) / (cm  * delta);
-				T dmdl = (fmu - fmu) / (cm * delta);// absurd if not spherical!
-				T dmdt = (fmv - fmv) / (cm * delta);
+				T dmdl = (fmup - fmu) / (cm * delta);// absurd if not spherical!
+				T dmdt = (fmvp - fmv) / (cm * delta);
 				T fG = vvi * dmdl - uui * dmdt;
 				XAdv.dhu[i] = (XFlux.Fqux[i] + XFlux.Fquy[i] - XFlux.Su[iright] - XFlux.Fquy[itop]) * cmdinv + fc * hi * vvi;
 				XAdv.dhv[i] = (XFlux.Fqvy[i] + XFlux.Fqvx[i] - XFlux.Sv[itop] - XFlux.Fqvx[iright]) * cmdinv - fc * hi * uui;

src/Advection.h

@@ -6,6 +6,7 @@
 #include "Arrays.h"
 #include "Forcing.h"
 #include "MemManagement.h"
+#include "Spherical.h"
 
 template <class T> __global__ void updateEVGPU(Param XParam, BlockP<T> XBlock, EvolvingP<T> XEv, FluxP<T> XFlux, AdvanceP<T> XAdv);
 template <class T> __host__ void updateEVCPU(Param XParam, BlockP<T> XBlock, EvolvingP<T> XEv, FluxP<T> XFlux, AdvanceP<T> XAdv);

src/ConserveElevation.cu

@@ -787,7 +787,7 @@ template <class T> void conserveElevationGHLeft(Param XParam, int ib, int ibLB,
 {
 	int ibn;
 	int ihalo, jhalo, ip, jp, iq, jq;
-	T delta = calcres(T(XParam.dx), XBlock.level[ib]);
+	T delta = calcres(T(XParam.delta), XBlock.level[ib]);
 	ihalo = -1;
 	ip = 0;
 
@@ -871,7 +871,7 @@ template <class T> __global__ void conserveElevationGHLeft(Param XParam, BlockP<
 	int ip, jp, iq, jq;
 
 	int ihalo, jhalo, ibn;
-	T delta = calcres(XParam.dx, lev);
+	T delta = calcres(XParam.delta, lev);
 
 
 	ihalo = -1;
@@ -926,7 +926,7 @@ template <class T> void conserveElevationGHRight(Param XParam, int ib, int ibRB,
 {
 	int ibn;
 	int ihalo, jhalo, ip, jp, iq, jq;
-	T delta = calcres(T(XParam.dx), XBlock.level[ib]);
+	T delta = calcres(T(XParam.delta), XBlock.level[ib]);
 	ihalo = XParam.blkwidth;
 	ip = XParam.blkwidth-1;
 
@@ -1004,7 +1004,7 @@ template <class T> __global__ void conserveElevationGHRight(Param XParam, BlockP
 
 	int ihalo, jhalo, iq, jq, ip, jp, ibn;
 
-	T delta = calcres(XParam.dx, lev);
+	T delta = calcres(XParam.delta, lev);
 
 	ihalo = blockDim.y;
 	jhalo = iy;
@@ -1061,7 +1061,7 @@ template <class T> void conserveElevationGHTop(Param XParam, int ib, int ibTL, i
 {
 	int ibn;
 	int ihalo, jhalo, ip, jp, iq, jq;
-	T delta = calcres(T(XParam.dx), XBlock.level[ib]);
+	T delta = calcres(T(XParam.delta), XBlock.level[ib]);
 	jhalo = XParam.blkwidth;
 	jp = XParam.blkwidth - 1;
 
@@ -1137,7 +1137,7 @@ template <class T> __global__ void conserveElevationGHTop(Param XParam, BlockP<T
 
 
 	int ihalo, jhalo, iq, jq, ip, jp, ibn;
-	T delta = calcres(XParam.dx, lev);
+	T delta = calcres(XParam.delta, lev);
 
 	ihalo = ix;
 	jhalo = iy+1;
@@ -1192,7 +1192,7 @@ template <class T> void conserveElevationGHBot(Param XParam, int ib, int ibBL, i
 {
 	int ibn;
 	int ihalo, jhalo, ip, jp, iq, jq;
-	T delta = calcres(T(XParam.dx), XBlock.level[ib]);
+	T delta = calcres(T(XParam.delta), XBlock.level[ib]);
 	jhalo = -1;
 	jp = 0;
 
@@ -1269,7 +1269,7 @@ template <class T> __global__ void conserveElevationGHBot(Param XParam, BlockP<T
 	int ip, jp, iq, jq;
 
 	int ihalo, jhalo, ibn;
-	T delta = calcres(XParam.dx, lev);
+	T delta = calcres(XParam.delta, lev);
 
 	ihalo = ix;
 	jhalo = -1;

src/FlowGPU.cu

@@ -34,15 +34,18 @@ template <class T> void FlowGPU(Param XParam, Loop<T>& XLoop, Forcing<float> XFo
 		cudaStreamDestroy(atmpstreams[0]);
 
 		//Calc dpdx and dpdy
-		gradient<< < gridDim, blockDim, 0 >> > (XParam.halowidth, XModel.blocks.active, XModel.blocks.level, (T)XParam.theta, (T)XParam.dx, XModel.Patm, XModel.datmpdx, XModel.datmpdy);
+
+		gradient << < gridDim, blockDim, 0 >> > (XParam.halowidth, XModel.blocks.active, XModel.blocks.level, (T)XParam.theta, (T)XParam.delta, XModel.Patm, XModel.datmpdx, XModel.datmpdy);
+
+
 		CUDA_CHECK(cudaDeviceSynchronize());
 		gradientHaloGPU(XParam, XModel.blocks, XModel.Patm, XModel.datmpdx, XModel.datmpdy);
 		//
 
 
 		refine_linearGPU(XParam, XModel.blocks, XModel.Patm, XModel.datmpdx, XModel.datmpdy);
 
-		gradient << < gridDim, blockDim, 0 >> > (XParam.halowidth, XModel.blocks.active, XModel.blocks.level, (T)XParam.theta, (T)XParam.dx, XModel.Patm, XModel.datmpdx, XModel.datmpdy);
+		gradient << < gridDim, blockDim, 0 >> > (XParam.halowidth, XModel.blocks.active, XModel.blocks.level, (T)XParam.theta, (T)XParam.delta, XModel.Patm, XModel.datmpdx, XModel.datmpdy);
 		CUDA_CHECK(cudaDeviceSynchronize());
 
 		gradientHaloGPU(XParam, XModel.blocks, XModel.Patm, XModel.datmpdx, XModel.datmpdy);