Fix Dynamic forcing

CyprienBosserelle · Jan 16, 2025 · bda4751 · bda4751
1 parent 306cc11
commit bda4751
Show file tree

Hide file tree

Showing 4 changed files with 14 additions and 8 deletions.
diff --git a/src/GridManip.cu b/src/GridManip.cu
@@ -453,7 +453,7 @@ template void InterpstepCPU<float, double>(int nx, int ny, int hdstep, double to
 template void InterpstepCPU<double, double>(int nx, int ny, int hdstep, double totaltime, double hddt, double*& Ux, double* Uo, double* Un);
 
 
-template <class T> __global__ void InterpstepGPU(int nx, int ny, int hdstp, T totaltime, T hddt, T*Ux, T* Uo, T* Un)
+template <class T> __global__ void InterpstepGPU(int nx, int ny, T totaltime, T beforetime, T aftertime, T*Ux, T* Uo, T* Un)
 {
 	unsigned int ix = blockIdx.x * blockDim.x + threadIdx.x;
 	unsigned int iy = blockIdx.y * blockDim.y + threadIdx.y;
@@ -464,19 +464,21 @@ template <class T> __global__ void InterpstepGPU(int nx, int ny, int hdstp, T to
 	__shared__ T Uxn[16][16];
 	//	__shared__ float Ums[16];
 
-
+	T hddt = aftertime - beforetime;
 
 	if (ix < nx && iy < ny)
 	{
 		Uxo[tx][ty] = Uo[ix + nx * iy]/**Ums[tx]*/;
 		Uxn[tx][ty] = Un[ix + nx * iy]/**Ums[tx]*/;
 
-		Ux[ix + nx * iy] = Uxo[tx][ty] + (totaltime - hddt * hdstp) * (Uxn[tx][ty] - Uxo[tx][ty]) / hddt;
+		Ux[ix + nx * iy] = Uxo[tx][ty] + (totaltime - beforetime) * (Uxn[tx][ty] - Uxo[tx][ty]) / (hddt);
 	}
 }
 //template __global__ void InterpstepGPU<int>(int nx, int ny, int hdstp, T totaltime, T hddt, T* Ux, T* Uo, T* Un);
-template __global__ void InterpstepGPU<float>(int nx, int ny, int hdstp, float totaltime, float hddt, float* Ux, float* Uo, float* Un);
-template __global__ void InterpstepGPU<double>(int nx, int ny, int hdstp, double totaltime, double hddt, double* Ux, double* Uo, double* Un);
+template __global__ void InterpstepGPU<float>(int nx, int ny, float totaltime, float beforetime, float aftertime, float* Ux, float* Uo, float* Un);
+template __global__ void InterpstepGPU<double>(int nx, int ny, double totaltime, double beforetime, double aftertime, double* Ux, double* Uo, double* Un);
+
+
 
 template <class T> void Copy2CartCPU(int nx, int ny, T* dest, T* src)
 {

diff --git a/src/GridManip.h b/src/GridManip.h
@@ -23,7 +23,8 @@ template <class T, class F> T interp2BUQ(T x, T y, F forcing);
 template <class T, class F> T interp2BUQ(T x, T y, T dx, F forcing);
 
 template <class T, class F> void InterpstepCPU(int nx, int ny, int hdstep, F totaltime, F hddt, T*& Ux, T* Uo, T* Un);
-template <class T> __global__ void InterpstepGPU(int nx, int ny, int hdstp, T totaltime, T hddt, T* Ux, T* Uo, T* Un);
+//template <class T> __global__ void InterpstepGPU(int nx, int ny, int hdstp, T totaltime, T hddt, T* Ux, T* Uo, T* Un);
+template <class T> __global__ void InterpstepGPU(int nx, int ny, T totaltime, T beforetime, T aftertime, T* Ux, T* Uo, T* Un);
 
 template <class T> void Copy2CartCPU(int nx, int ny, T* dest, T* src);
 

diff --git a/src/Read_netcdf.cu b/src/Read_netcdf.cu
@@ -546,6 +546,7 @@ int readnctime2(int ncid,char * timecoordname,std::string refdate,size_t nt, dou
 	for (int it = 0; it < nt; it++)
 	{
 		time[it] = time[it] * fac + offset;
+		//printf("%f\n", time[it]);
 	}
 
 	///status = nc_close(ncid);

diff --git a/src/Updateforcing.cu b/src/Updateforcing.cu
@@ -98,7 +98,9 @@ void Forcingthisstep(Param XParam, double totaltime, DynForcingP<float> &XDynFor
 		// Interpolate the forcing array to this time 
 		if (XParam.GPUDEVICE >= 0)
 		{
-			InterpstepGPU << <gridDimDF, blockDimDF, 0 >> > (XDynForcing.nx, XDynForcing.ny, XDynForcing.instep - 1, float(totaltime), float(XDynForcing.dt), XDynForcing.now_g, XDynForcing.before_g, XDynForcing.after_g);
+			float bftime = float(XDynForcing.to+XDynForcing.dt*(XDynForcing.instep-1));
+			float aftime = float(XDynForcing.to + XDynForcing.dt * (XDynForcing.instep));
+			InterpstepGPU << <gridDimDF, blockDimDF, 0 >> > (XDynForcing.nx, XDynForcing.ny, float(totaltime), bftime,aftime, XDynForcing.now_g, XDynForcing.before_g, XDynForcing.after_g);
 			CUDA_CHECK(cudaDeviceSynchronize());
 
 			CUDA_CHECK(cudaMemcpyToArray(XDynForcing.GPU.CudArr, 0, 0, XDynForcing.now_g, XDynForcing.nx * XDynForcing.ny * sizeof(float), cudaMemcpyDeviceToDevice));
@@ -377,7 +379,7 @@ template <class T> __global__ void AddrainforcingImplicitGPU(Param XParam, Loop<
 
 
 	Rainhh = max(Rainhh / T(1000.0) / T(3600.0) * T(XLoop.dt), T(0.0)) * XBlock.activeCell[i]; // convert from mm/hrs to m/s and 
-
+	//printf("%f\n", Rainhh);
 	T qvol = hi / (hi + Rainhh);
 
 	XEv.h[i] = hi + Rainhh;