From c0d4bdedbe4984571219b2efd3dbb8e98fe80431 Mon Sep 17 00:00:00 2001 From: leogean3 <61390859+leogean3@users.noreply.github.com> Date: Sun, 2 May 2021 16:42:25 -0600 Subject: [PATCH] Create WDPM_mulGPU.c Upload the multiple-GPU code of WDPM --- WDPM_mulGPU.c | 2440 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 2440 insertions(+) create mode 100644 WDPM_mulGPU.c diff --git a/WDPM_mulGPU.c b/WDPM_mulGPU.c new file mode 100644 index 0000000..98ccb0e --- /dev/null +++ b/WDPM_mulGPU.c @@ -0,0 +1,2440 @@ +#define PROGRAM_FILE "runoff.cl" +#define KERNEL_FUNC1 "add" +#define KERNEL_FUNC2 "subtract" +#define KERNEL_FUNC3 "ddrain" + +#define CL_USE_DEPRECATED_OPENCL_1_1_APIS + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#define max( a, b ) ( ((a) > (b)) ? (a) : (b) ) +#define min( a, b ) ( ((a) < (b)) ? (a) : (b) ) + +#ifdef __APPLE__ +#include +#else +#include "CL/cl.h" +#endif + + + +#ifdef _WIN32 +#include +#include /* obtain _kbhit() and _getch() for killswitch - Ryan Spies */ + +#if defined(_MSC_VER) || defined(_MSC_EXTENSIONS) +#define DELTA_EPOCH_IN_MICROSECS 11644473600000000Ui64 +#else +#define DELTA_EPOCH_IN_MICROSECS 11644473600000000ULL +#endif + +struct timeZone +{ + int tz_minuteswest; /* minutes W of Greenwich */ + int tz_dsttime; /* type of dst correction */ +}; + +int gettimeofday(struct timeval *tv, struct timeZone *tz) +{ + FILETIME ft; + unsigned __int64 tmpres = 0; + static int tzflag; + + if (NULL != tv) + { + GetSystemTimeAsFileTime(&ft); + + tmpres |= ft.dwHighDateTime; + tmpres <<= 32; + tmpres |= ft.dwLowDateTime; + + /*converting file time to unix epoch*/ + tmpres -= DELTA_EPOCH_IN_MICROSECS; + tmpres /= 10; /*convert into microseconds*/ + tv->tv_sec = (long)(tmpres / 1000000UL); + tv->tv_usec = (long)(tmpres % 1000000UL); + } + + if (NULL != tz) + { + if (!tzflag) + { + _tzset(); + tzflag++; + } + tz->tz_minuteswest = _timezone / 60; + tz->tz_dsttime = _daylight; + } + return 0; +} +#else +// If we aren't using Windows, get getch() and kbhit() this way - Ryan Spies +#include +#include +#include + +struct termios orig_termios; + +void reset_terminal_mode() +{ + tcsetattr(0, TCSANOW, &orig_termios); +} + +void set_conio_terminal_mode() +{ + struct termios new_termios; + + /* take two copies - one for now, one for later */ + tcgetattr(0, &orig_termios); + memcpy(&new_termios, &orig_termios, sizeof(new_termios)); + + /* register cleanup handler, and set the new terminal mode */ + atexit(reset_terminal_mode); + cfmakeraw(&new_termios); + tcsetattr(0, TCSANOW, &new_termios); +} + +int kbhit() +{ + struct timeval tv = { 0L, 0L }; + fd_set fds; + FD_ZERO(&fds); + FD_SET(0, &fds); + return select(1, &fds, NULL, NULL, &tv); +} + +int getch() +{ + int r; + unsigned char c; + if ((r = read(0, &c, sizeof(c))) < 0) { + return r; + } else { + return c; + } +} +#endif + +cl_uint get_num_devices(int opencl) { + cl_platform_id *platforms; + int err; + cl_uint num_platforms, num_devices; //must be uint + + /* Identify all platforms platforms */ + clGetPlatformIDs(0, NULL, &num_platforms); + // printf("There are %d platforms \n", num_platforms); + platforms = (cl_platform_id*) malloc (num_platforms*sizeof(cl_platform_id)); + err = clGetPlatformIDs(2, platforms, &num_platforms); + + if(err < 0) { + perror("Couldn't identify a platform"); + exit(1); + } + + if (opencl == 0){ + err = clGetDeviceIDs(platforms[opencl], CL_DEVICE_TYPE_GPU, 0, NULL, &num_devices); + if(err == CL_DEVICE_NOT_FOUND) { + err = clGetDeviceIDs(platforms[0], CL_DEVICE_TYPE_GPU, 0, NULL, &num_devices); + } + if(err < 0) { + perror("Couldn't access any devices"); + exit(1); + } +} + + if (opencl == 1){ + err = clGetDeviceIDs(platforms[0], CL_DEVICE_TYPE_CPU, 0, NULL, &num_devices); + if(err == CL_DEVICE_NOT_FOUND) { + err = clGetDeviceIDs(platforms[0], CL_DEVICE_TYPE_CPU, 0, NULL, &num_devices); + } + if(err < 0) { + perror("Couldn't access any devices"); + exit(1); + } + } + + return num_devices; +} + + /* Find a GPU or CPU associated with the first available platform */ +cl_device_id* create_devices(int opencl, cl_uint num_devices) { + cl_platform_id *platforms; + cl_device_id* dev; + int err; + cl_uint num_platforms; //must be uint + + /* Identify all platforms platforms */ + clGetPlatformIDs(0, NULL, &num_platforms); + // printf("There are %d platforms \n", num_platforms); + platforms = (cl_platform_id*) malloc (num_platforms*sizeof(cl_platform_id)); + err = clGetPlatformIDs(2, platforms, &num_platforms); + dev = (cl_device_id *)malloc(num_devices * sizeof(cl_device_id) ); + + if(err < 0) { + perror("Couldn't identify a platform"); + exit(1); + } + + // Choose all available devices of one type - Ryan Spies + if (opencl == 0){ + err = clGetDeviceIDs(platforms[opencl], CL_DEVICE_TYPE_GPU, num_devices, dev, NULL); + if(err == CL_DEVICE_NOT_FOUND) { + err = clGetDeviceIDs(platforms[0], CL_DEVICE_TYPE_GPU, num_devices, dev, NULL); + } + if(err < 0) { + perror("Couldn't access any devices"); + exit(1); + } + } + + if (opencl == 1){ + err = clGetDeviceIDs(platforms[0], CL_DEVICE_TYPE_CPU, 1, dev, NULL); + if(err == CL_DEVICE_NOT_FOUND) { + err = clGetDeviceIDs(platforms[opencl-1], CL_DEVICE_TYPE_CPU, 1, dev, NULL); + } + if(err < 0) { + perror("Couldn't access any devices"); + exit(1); + } + } + + return dev; +} + +/* Create program from a file and compile it */ +cl_program build_program(cl_context ctx, cl_device_id* dev, cl_uint num_devices, const char* filename) { + + cl_program program; + FILE *program_handle; + char *program_buffer, *program_log; + size_t program_size, log_size; + int err, i; + + /* Read program file and place content into buffer */ + program_handle = fopen(filename, "r"); + if(program_handle == NULL) { + perror("Couldn't find the program file"); + exit(1); + } + fseek(program_handle, 0, SEEK_END); + program_size = ftell(program_handle); + rewind(program_handle); + program_buffer = (char*)malloc(program_size + 1); + program_buffer[program_size] = '\0'; + fread(program_buffer, sizeof(char), program_size, program_handle); + fclose(program_handle); + + /* Create program from file */ + program = clCreateProgramWithSource(ctx, 1, (const char**)&program_buffer, &program_size, &err); + if(err < 0) { + perror("Couldn't create the program"); + exit(1); + } + free(program_buffer); + + /* Build program */ + err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL); + if(err < 0) { + /* Find size of log and print to std output */ + for(i=0;i 1) + { + for(i=0;i missingvalue ){ + basincount++; + } + } + } + char *temp1, *temp2; + if (strcmp(activity, "add") == 0) { + // get water volume at beginning + initial_vol=0; + for(i=0; imissingvalue){ + initial_vol+=water[i][j]; + } + } + } + initial_vol=initial_vol*cellarea; + + temp1 = upcase(ScratchFileName); + // check for existence of scratch file & read if available + if ( strcmp(temp1,"NULL")!=0 ){ + if ( file_exist(ScratchFileName)){ + printf("%s\n"," "); + printf("%30s\n", "Scratch file found"); + read_water_array(ScratchFileName, numrows, numcols); + } + else + { + printf("%s\n"," "); + printf("%30s\n", "No Scratch file found"); + printf("%s\n"," "); + printf("%30s\n", "New Scratch will be saved"); + printf("%s\n"," "); + printf("%30s\n", "Now proceeding with Waterfile checking"); + printf("%s\n", " "); + // check for water values + temp2 = upcase(WaterFileName); + if (strcmp(temp2,"NULL") != 0){ + if (file_exist(WaterFileName)){ + // file exists, so read values + printf("%30s\n", "Existing water file found"); + read_water_array(WaterFileName, numrows, numcols); + // get water volume at beginning + initial_vol = 0; + for (i = 0; imissingvalue){ + initial_vol += water[i][j]; + } + } + } + initial_vol = initial_vol*cellarea; + } + else + { + // no water file + printf("%30s\n", "Water file missing, will be created"); + for(i=0; imissingvalue){ + water[i][j]=0; + } + } + } + } + } + else + { + printf("%30s\n", "Water file will be created"); + for(i=0; imissingvalue){ + water[i][j]=0; + } + } + } + } + free(temp2); + // add water to array + for(i=0; imissingvalue && water[i][j]>0){ + water[i][j]+=addwater; + } + } + } + for(i=0; imissingvalue && water[i][j]<=0){ + water[i][j]=addwater*rof; + } + } + } + } + } + else + { + // check for water values + temp2 = upcase(WaterFileName); + if (strcmp(temp2,"NULL") != 0){ + if (file_exist(WaterFileName)){ + // file exists, so read values + printf("%30s\n", "Existing water file found"); + read_water_array(WaterFileName, numrows, numcols); + } + else + { + // no water file + printf("%30s\n", "Water file missing, will be created"); + for(i=0; imissingvalue){ + water[i][j]=0; + } + } + } + } + } + else + { + printf("%30s\n", "Water file will be created"); + for(i=0; imissingvalue){ + water[i][j]=0; + } + } + } + } + free(temp2); + // add water to array + for(i=0; imissingvalue && water[i][j]>0){ + water[i][j]+=addwater; + } + } + } + for(i=0; imissingvalue && water[i][j]<=0){ + water[i][j]=addwater*rof; + } + } + } + } + free(temp1); + // intialize big arrays + for(i=0; i0){ + initial_vol += water[i][j]; + } + } + } + initial_vol = initial_vol*cellarea; + temp1 = upcase(ScratchFileName); + // check for existence of scratch file & read if available + if ( strcmp(temp1,"NULL")!=0 ){ + if ( file_exist(ScratchFileName)){ + printf("%s\n"," "); + printf("%30s\n", "Scratch file found"); + read_water_array(ScratchFileName, numrows, numcols); + } + else + { + printf("%s\n"," "); + printf("%30s\n", "No Scratch file found"); + printf("%s\n"," "); + printf("%30s\n", "New Scratch will be saved."); + printf("%s\n"," "); + printf("%30s\n", "Now proceeding with Waterfile checking"); + printf("%s\n", " "); + temp2 = upcase(WaterFileName); + // check for water values + if (strcmp(temp2,"NULL") != 0){ + if ( file_exist(WaterFileName)){ + // file exists, so read values + printf("%30s\n", "Existing water file found"); + read_water_array(WaterFileName, numrows, numcols); + // get water volume at beginning + initial_vol = 0; + for (i = 0; i0){ + initial_vol += water[i][j]; + } + } + } + initial_vol = initial_vol*cellarea; + } + else + { + // no water file + printf("%30s\n", "Water file missing, will be created"); + for(i=0; imissingvalue){ + water[i][j]=max(water[i][j]-subtractwater,0); + } + } + } + } + } + else + { + temp2 = upcase(WaterFileName); + // check for water values + if (strcmp(temp2,"NULL") != 0){ + if ( file_exist(WaterFileName)){ + // file exists, so read values + printf("%30s\n", "Existing water file found"); + read_water_array(WaterFileName, numrows, numcols); + } + else + { + // no water file + printf("%30s\n", "Water file missing, will be created"); + for(i=0; imissingvalue){ + water[i][j]=max(water[i][j]-subtractwater,0); + } + } + } + } + free(temp1); + // intialize big arrays + for(i=0; i0){ + if(bigdem[i][j] 0) + { + for(int i=0; i(offsets[i]+i)) && (draincolmissingvalue){ + initial_vol+=water[i][j]; + } + } + } + initial_vol=initial_vol*cellarea; + totaldrain=max(bigwater[drainrow][draincol],0); + basin_area = basincount*cellarea; + + print_basin_summary(basin_area, initial_vol, drainrow, draincol, bigdem[drainrow][draincol]); + print_iteration_summary_headings(); + } + + + + // now do water runoff until finished + done = false; + k = 0; + int oi, oj; + // do a set of iterations and then test for convergence + gettimeofday(&starttime, NULL); + while (done == false){ + for(i=0; i<(numrows+2); i++) + { + for(j=0; j<(numcols+2); j++) + { + if(bigwater[i][j] < thres) + { + bigwater[i][j] = 0; + } + } + } + if (strcmp(activity, "drain") == 0) { + olddrain = totaldrain; + } + for(i=0; i 0.0 && (bigdem[row][col] > missingvalue) + && (row != drainrow || col !=draincol )){ + runoffd(row, col, drainrow, draincol, missingvalue); + } + } + } + } + } + totaldrain = totaldrain + drain(drainrow,draincol); + }// 1000 iterations + } + else if (strcmp(activity, "add") == 0) + { + for (i=0; i 0.0 && (bigdem[row][col] > missingvalue)){ + runoffs(row, col, missingvalue); + } + } + } + } + } + }// 1000 iterations + } + else if (strcmp(activity, "subtract") == 0) + { + for (i=0; i 0.0 && (bigdem[row][col] > missingvalue)){ + runoffs(row, col, missingvalue); + } + } + } + } + } + }// 1000 iterations + } + k = k + IterationNum; + } + else if(cpu == 1) + { + + if(num_devices > 1) + { + // Flattened the Matrix to 1D array + for(i=0; i 1) + { + if(i>0) + { + for(j=0; j<(num_devices-1); j++) + { + err = clEnqueueWriteBuffer(queues[j], d_water[j], CL_FALSE, (section_sizes[j]-2)*(numrows+2)*sizeof(double), (numrows+2)*2*sizeof(double), b_dem[j], 0, NULL, &event); + exitOnFail(err, "load updated subwater buffer"); + err = clWaitForEvents(1, &event); + exitOnFail(err, "wait for updated subwater buffer to load"); + clReleaseEvent(event); + } + for(j=1; j0.0 && (bigdem[m][n]>missingvalue)) + { + runoffa(m, n, missingvalue); + } + } + } + } + for(int o=0; o<(num_devices-1); o++) + { + for(int m=0; m<2; m++) + { + for(int n=0; n<(numrows+2); n++) + { + b_dem[o][m*(numrows+2)+n]=bigwater[n][offsets[o+1]-1+o+m]; + b_water[o+1][m*(numrows+2)+n]=bigwater[n][offsets[o+1]-1+o+m+2]; + } + } + } + + if (strcmp(activity, "drain") == 0) { + totaldrain = 0; + int count; + for (count=0; countmissingvalue){ + if(diff[i][j]>max_diff){ + max_diff = diff[i][j]; + } + } + } + } + + if (strcmp(activity, "drain") == 0) { + diffdrain = (fabs(totaldrain - olddrain))*cellarea; + final_vol=0; + for(i=0; imissingvalue){ + final_vol+=bigwater[i][j]; + } + } + } + final_vol=final_vol*cellarea; + } + gettimeofday(¤ttime, NULL); + if (strcmp(activity, "drain") == 0) { + printf("%7s %d %7s %8.6f %5s %10.1f %5s %12.1f %5s %8.2f\n", "", k, "", max_diff,"", diffdrain, "", final_vol, "", + (double)(currenttime.tv_usec-starttime.tv_usec)/1000000 + (double)(currenttime.tv_sec-starttime.tv_sec)); + } + else + { + printf("%7s %d %7s %8.6f %5s %8.2f\n", "", k, "", + max_diff,"", (double)(currenttime.tv_usec-starttime.tv_usec)/1000000 + (double)(currenttime.tv_sec-starttime.tv_sec)); + } + + // check to see if graceful exit triggered - Ryan Spies + #ifdef _WIN32 + if (_kbhit()) + { + exitChar = _getch(); + if (exitChar == 'q') + { + gracefulExit = true; + printf("Now beginning graceful exit...\n"); + } + } + #else + if (kbhit()) + { + exitChar = getch(); + if (exitChar == 'q') + { + gracefulExit = true; + printf("Now beginning graceful exit...\n"); + } + } + #endif + // check to see if finished + if (strcmp(activity, "drain") == 0) { + temp1 = upcase(ScratchFileName); + if (max_diff <= elevation_tolerance || diffdrain < drain_tolerance || ((k == iteration_limit) && isKillSwitch) || gracefulExit) { + done = true; + } + else if (strcmp(temp1, "NULL")!=0) { + // write output to ArcGIS file + // convert back to small arrays + for(i=0; i0.001 && dem[i][j]> missingvalue){ + watercount++; + } + } + } + + watertotal=0; + for(i=0; i missingvalue){ + watertotal+=water[i][j]; + } + } + } + + final_vol=0; + for(i=0; imissingvalue){ + final_vol+=water[i][j]; + } + } + } + final_vol=final_vol*cellarea; + meanwater = watertotal/((float)watercount); + basincount=0; + for(i=0; imissingvalue){ + basincount++; + } + } + } + + waterfrac = (float)watercount/(float)basincount; + + if (strcmp(activity, "drain") == 0) { + drainvol = totaldrain*cellarea; + draindepth = (drainvol/((float)basincount*cellarea))*1000; + } + + maxdepth=water[0][0]; + for(i=0; imaxdepth){ + maxdepth = water[i][j]; + } + } + } + maxdepth = maxdepth*1000; + + if (strcmp(activity, "drain") == 0) { + print_results(initial_vol, final_vol, drainvol, waterfrac, meanwater, draindepth, maxdepth, iteration_limit); + } + else + { + print_results(initial_vol, final_vol, 0, waterfrac, meanwater, 0, maxdepth, iteration_limit); + } + + //write output to ArcGIS file + write_gis(OutputFileName, numrows, numcols, dem_header_line, dem_header_value); + + // print execution time + gettimeofday(¤ttime, NULL); + printf("%20s %10.2f %s\n","Run Time", (double)(currenttime.tv_usec-starttime.tv_usec)/1000000 + (double)(currenttime.tv_sec-starttime.tv_sec), "s"); + if(cpu == 1) + { + for(i=0;i 1){ + free(b_water[i]); + free(b_dem[i]); + } + } + err= clReleaseProgram(program); + err= clReleaseContext(context); + free(bbigdem); + free(bbigwater); + } + for (i=0; i<7; i++){ + free(dem_header_line[i]); + } + free(dem_header_line); + for (i=0; i=97 && string[i]<=122){ + USTR[i]=string[i]-32; + } + else + { + USTR[i]=string[i]; + } + } + return USTR; +} + + + +void write_gis(char filename[80], int numrows, int numcols, + char* header_line[8], double header_value[8]){ + //write output to ArcGIS file + FILE *fgis; + fgis=fopen(filename,"w"); + fprintf(fgis, "%s %d\n", header_line[1], (int)header_value[1]); + fprintf(fgis, "%s %d\n", header_line[2], (int)header_value[2]); + fprintf(fgis, "%s %14.6f\n", header_line[3], header_value[3]); + fprintf(fgis, "%s %14.6f\n", header_line[4], header_value[4]); + fprintf(fgis, "%s %9.6f\n", header_line[5], header_value[5]); + fprintf(fgis, "%s %14.6f\n", header_line[6], header_value[6]); + int row=0; + int col=0; + for (row=0; row < numrows; row++){ + for (col=0; col < numcols; col++) + { + fprintf(fgis, "%f ", water[row][col]); + } + fprintf(fgis, "\n"); + } + fclose(fgis); +} + +void read_water_array(char filename[80], int numrows, int numcols) { + // reads data from ArcGIS file + // 1st read and discard header + FILE *fwarray; + fwarray=fopen(filename,"r"); + int row; + int col; + fscanf(fwarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fwarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fwarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fwarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fwarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fwarray,"%s %lf", dem_header_name, &dem_header_value[0]); + for (row=0; row < numrows; row++){ + for (col=0; col < numcols; col++) + { + fscanf(fwarray,"%lf",&water[row][col]); + } + } + fclose(fwarray); +} + + +void read_dem_array(char filename[80], int numrows, int numcols) { + // reads data from ArcGIS file + // 1st read and discard header + FILE *fdarray; + fdarray=fopen(filename,"r"); + int row; + int col; + fscanf(fdarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fdarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fdarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fdarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fdarray,"%s %lf", dem_header_name, &dem_header_value[0]); + fscanf(fdarray,"%s %lf", dem_header_name, &dem_header_value[0]); + for (row=0; row < numrows; row++){ + for (col=0; col < numcols; col++) + { + fscanf(fdarray,"%lf",&dem[row][col]); + } + } + fclose(fdarray); +} + + +void read_gis_header(char filename[80], bool print_output, char* header_line[8], double header_value[8]){ + + if (print_output){ + printf("%s\n", " "); + printf("%30s\n", "ArcGIS file header"); + printf("%30s %d\n", header_line[1], (int)header_value[1]); + printf("%30s %d\n", header_line[2], (int)header_value[2]); + printf("%30s %9.1f\n", header_line[3], header_value[3]); + printf("%30s %9.1f\n", header_line[4], header_value[4]); + printf("%30s %9.1f\n", header_line[5], header_value[5]); + printf("%30s %9.1f\n", header_line[6], header_value[6]); + } +} + + +void print_copyright(char *activity){ + printf("%s\n", " "); + printf("%s\n", " "); + if (strcmp(activity, "add") == 0) { + printf("%s\n", "WDPM_add 1.0 - Wetland DEM Ponding Model - parallel version"); + } + else if (strcmp(activity, "subtract") == 0) { + printf("%s\n", "WDPM_subtract 1.0 - Wetland DEM Ponding Model - parallel version"); + } + else if (strcmp(activity, "drain") == 0) { + printf("%s\n", "WDPM_drain 1.0 - Wetland DEM Ponding Model - parallel version"); + } + printf("%s\n", "Copyright (C) 2010,2012 Kevin Shook, Centre for Hydrology"); + printf("%s\n", "--------------------------------------------------------------------"); + printf("%s\n", " "); + printf("%s\n", "This program is free software: you can redistribute it and/or modify"); + printf("%s\n", "it under the terms of the GNU General Public License as published by"); + printf("%s\n", "the Free Software Foundation, either version 3 of the License, or"); + printf("%s\n", "(at your option) any later version."); + printf("%s\n", " "); + printf("%s\n", "This program is distributed in the hope that it will be useful,"); + printf("%s\n", "but WITHOUT ANY WARRANTY; without even the implied warranty of"); + printf("%s\n", "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the"); + printf("%s\n", "GNU General Public License for more details."); + printf("%s\n", " "); + printf("%s\n", "You should have received a copy of the GNU General Public License"); + printf("%s\n", "along with this program. If not, see ."); + printf("%s\n", " "); + + if (strcmp(activity, "add") == 0) { + printf("%s\n", "This program adds water to an ArcGIS ASCII file of water runoff"); + printf("%s\n", "and redistributes water over the DEM"); + } + else if (strcmp(activity, "subtract") == 0) { + printf("%s\n", "This program removes a depth water to an ArcGIS ASCII file of water depths"); + printf("%s\n", "and redistributes water over the DEM"); + } + else if (strcmp(activity, "drain") == 0) { + printf("%s\n", "This program drains an ArcGIS ASCII file of water runoff"); + printf("%s\n", "from the lowest point in the DEM, which acts as a drain"); + } + printf("%s\n", "From the algorithm of Shapiro, M., & Westervelt, J. (1992). "); + printf("%s\n", "An Algebra for GIS and Image Processing (pp. 1-22)."); + printf("%s\n", " "); + printf("%s\n", " "); +} + + +void print_arg_list(){ + printf("%s\n", " "); + printf("%s\n", "Program arguments in order of specification"); + if (strcmp(activity,"add")==0){ + printf("%s\n", "Add module specified"); + } + else if (strcmp(activity,"subtract")==0){ + printf("%s\n", "Subtract module specified"); + } + else if (strcmp(activity,"drain")==0){ + printf("%s\n", "Drain module specified"); + } + printf("%s\n", "Path and Name of Report file"); + printf("%s\n", "DEM file name (string)"); + printf("%s\n", "Water file name (string)"); + printf("%s\n", "Output file name (string)"); + printf("%s\n", "Scratch file name (string) - Optional, use --NULL-- to omit"); + if (strcmp(activity,"add")==0){ + printf("%s\n", "Depth of water to add (mm) (real)"); + printf("%s\n", "Water runoff fraction (real)"); + printf("%s\n", "Elevation tolerance (mm) (real)"); + } + else if (strcmp(activity,"subtract")==0){ + printf("%s\n", "Depth of water to remove (mm) (real)"); + printf("%s\n", "Elevation tolerance (mm) (real)"); + } + else if (strcmp(activity,"drain")==0){ + printf("%s\n", "Elevation tolerance (mm) (real)"); + printf("%s\n", "Drain tolerance (m3) (real)"); + } + printf("%s\n", "Specify 0 for serial CPU and 1 for opencl "); + printf("%s\n", "Specify 0 for OpenCL CPU and 1 for opencl GPU "); + printf("%s\n", "Zero depth threshold (mm) (real)"); + printf("%s\n", "Iteration limit (multiple of 1000)"); + printf("%s\n", " "); +} + +void print_narg_list(){ + printf("%s\n", "Module name: add"); + printf("%s\n", "Path and Name of Report file"); + printf("%s\n", "DEM file name (string)"); + printf("%s\n", "Water file name (string)"); + printf("%s\n", "Output file name (string)"); + printf("%s\n", "Scratch file name (string) - Optional, use --NULL-- to omit"); + printf("%s\n", "Depth of water to add (mm) (real)"); + printf("%s\n", "Water runoff fraction (real)"); + printf("%s\n", "Elevation tolerance (mm) (real)"); + printf("%s\n", "Specify 0 for serial CPU and 1 for opencl "); + printf("%s\n", "Specify 0 for OpenCL CPU and 1 for opencl GPU "); + printf("%s\n", "Zero depth threshold (mm) (real) "); + printf("%s\n", "Iteration limit (multiple of 1000)"); + printf("%s\n", " "); + printf("%s\n", " "); + printf("%s\n", "Module name: subtract"); + printf("%s\n", "Path and Name of Report file"); + printf("%s\n", "DEM file name (string)"); + printf("%s\n", "Water file name (string)"); + printf("%s\n", "Output file name (string)"); + printf("%s\n", "Scratch file name (string) - Optional, use --NULL-- to omit"); + printf("%s\n", "Depth of water to remove (mm) (real)"); + printf("%s\n", "Elevation tolerance (mm) (real)"); + printf("%s\n", "Specify 0 for serial CPU and 1 for opencl "); + printf("%s\n", "Specify 0 for OpenCL CPU and 1 for opencl GPU "); + printf("%s\n", "Zero depth threshold (mm) (real) "); + printf("%s\n", "Iteration limit (multiple of 1000)"); + printf("%s\n", " "); + printf("%s\n", " "); + printf("%s\n", "Module name: drain"); + printf("%s\n", "Path and Name of Report file"); + printf("%s\n", "DEM file name (string)"); + printf("%s\n", "Water file name (string)"); + printf("%s\n", "Output file name (string)"); + printf("%s\n", "Scratch file name (string) - Optional, use --NULL-- to omit"); + printf("%s\n", "Elevation tolerance (mm) (real)"); + printf("%s\n", "Drain tolerance (m3) (real)"); + printf("%s\n", "Specify 0 for serial CPU and 1 for opencl "); + printf("%s\n", "Specify 0 for OpenCL CPU and 1 for opencl GPU "); + printf("%s\n", "Zero depth threshold (mm) (real) "); + printf("%s\n", "Iteration limit (multiple of 1000)"); + printf("%s\n", " "); +} + + +void print_args (char demfile[80], char waterfile[80], char outputfile[80], + char scratchfile[80], double addwater, double subtractwater, double rof, + double eltol, double draintol, char *activity, int cpu, double thres, int iterlim){ + printf("%30s\n", "WDPM Parameters"); + printf("%30s %s\n", "Function used:", activity); + printf("%30s %s\n", "DEM file:", demfile); + printf("%30s %s\n", "Water file:", waterfile); + printf("%30s %s\n", "Output file:", outputfile); + printf("%30s %s\n", "Scratch file:", scratchfile); + if (strcmp(activity,"add")==0){ + printf("%30s %0.1f %s\n", "Water added:", addwater, "mm"); + printf("%30s %0.1f\n", "Runoff fraction:", rof); + printf("%30s %0.1f %s\n", "Elevation tolerance:", eltol, "mm"); + printf("%30s %0.4f %s\n", "Zero depth threshold:", thres, "mm"); + } + if (strcmp(activity,"subtract")==0){ + printf("%30s %0.1f %s\n", "Water subtracted:", subtractwater, "mm"); + printf("%30s %0.1f %s\n", "Elevation tolerance:", eltol, "mm"); + printf("%30s %0.4f %s\n", "Zero depth threshold:", thres, "mm"); + } + if (strcmp(activity,"drain")==0){ + printf("%30s %0.1f %s\n", "Elevation tolerance:", eltol, "mm"); + printf("%30s %0.1f %s\n", "Drain tolerance:", draintol, "m3"); + printf("%30s %0.4f %s\n", "Zero depth threshold:", thres, "mm"); + } + if (cpu==0){ + printf("%s\n", " "); + printf("%41s\n", "Using Serial CPU for Computation"); + } + else + { + printf("%s\n", " "); + printf("%41s\n", "Using Parallel OpenCL for Computation"); + } + if (iterlim == 0){ + printf("%41s\n", "No iteration limit specified"); + } + else + { + printf("%30s %i\n", "Iteration limit:", iterlim); + } +} + +void print_iteration_summary_headings(){ + if (strcmp(activity, "add") == 0) { + printf("%s\n", " "); + printf("%30s\n", "Doing calculations"); + printf("%15s %15s %15s\n", "iterations", "max diff", "run time"); + printf("%13s %14s %15s\n", " ", "(m)", "(s)"); + } + else if (strcmp(activity, "subtract") == 0) { + printf("%s\n", " "); + printf("%30s\n", "Doing calculations"); + printf("%15s %15s %15s\n", "iterations", "max diff", "run time"); + printf("%13s %14s %15s\n", " ", "(m)", "(s)"); + } + else if (strcmp(activity, "drain") == 0) { + printf("%s\n", " "); + printf("%30s\n", "Doing calculations"); + printf("%15s %15s %15s %15s %15s\n", "iterations", "max diff", "vol change", "water left", "run time"); + printf("%13s %14s %15s %16s %17s\n", " ", "(m)", "(m3)", "(m3)", "(s)"); + } +} + +void print_basin_summary(double basin_area, double initial_vol, int drainrow, int draincol, double minel){ + printf("%s\n", " "); + printf("%30s\n", "Basin summary"); + printf("%20s %10.1f %s\n", "Basin area:", basin_area, "m2"); + printf("%20s %10.1f %s\n", "Initial volume:", initial_vol, "m3"); + printf("%20s %d\n", "Drain column:", draincol); + printf("%20s %d\n", "Drain row:", drainrow); + printf("%20s %10.3f %s\n", "Min DEM elevation:", minel, "m"); +} + + +void print_results(double initial_vol, double final_vol, double drainvol, double waterfrac, + double meanwater, double draindepth, double maxdepth, int iterlim){ + if (strcmp(activity, "drain") == 0) { + printf("%s\n", " "); + printf("%30s\n", "WDPM run summary"); + printf("%20s %10.1f %s\n", "Initial volume",initial_vol,"m3"); + printf("%20s %10.1f %s\n", "Final volume",final_vol,"m3"); + printf("%20s %10.1f %s\n", "Volume change",initial_vol-final_vol,"m3"); + printf("%20s %10.1f %s\n", "Volume drained",drainvol,"m3"); + printf("%20s %10.3f %s\n", "Final water coverage", waterfrac,""); + printf("%20s %10.1f %s\n", "Mean water depth", meanwater*1000.,"mm"); + printf("%20s %10.1f %s\n", "Depth drained", draindepth,"mm "); + printf("%20s %10.1f %s\n", "Max water depth", maxdepth,"mm "); + printf("%20s %d\n", "Iteration limit", iterlim); + } + else + { + printf("%s\n", " "); + printf("%30s\n", "WDPM run summary"); + printf("%20s %10.1f %s\n", "Initial volume",initial_vol,"m3"); + printf("%20s %10.1f %s\n", "Final volume",final_vol,"m3"); + printf("%20s %10.1f %s\n", "Volume change",final_vol - initial_vol,"m3"); + printf("%20s %10.3f %s\n", "Final water coverage", waterfrac,""); + printf("%20s %10.1f %s\n", "Mean water depth", meanwater*1000.,"mm"); + printf("%20s %10.1f %s\n", "Max water depth", maxdepth,"mm "); + printf("%20s %d\n", "Iteration limit", iterlim); + } + +} + +double drain(int drainrow, int draincol) { + // drains water to specified location + int i,j; + double drainn; + double **waterslice, **demslice; + waterslice = malloc(3 * sizeof(double *)); + demslice = malloc(3 * sizeof(double *)); + for (i=0; i<3; i++){ + waterslice[i] = malloc(3 * sizeof(double)); + demslice[i] = malloc(3 * sizeof(double)); + } + // remove existing water + for (i=0; i<3; i++){ + for (j=0; j<3; j++){ + demslice[i][j] = bigdem[drainrow-1+i][draincol-1+j]; + waterslice[i][j] = bigwater[drainrow-1+i][draincol-1+j]; + } + } + drainn=0; + for (i=0; i<3; i++){ + for (j=0; j<3; j++){ + if(demslice[i][j] > missingvalue && waterslice[i][j] > 0){ + drainn+=waterslice[i][j]; + } + } + } + for (i=drainrow-1; i<=drainrow+1; i++){ + for (j=draincol-1; j<=draincol+1; j++){ + bigwater[i][j] = 0.0 ; + } + } + for (i=0; i<3; i++){ + free(waterslice[i]); + free(demslice[i]); + } + free(waterslice); + free(demslice); + return drainn; +} + + + +void runoffa(int centerrow, int centercol, double missingvalue) { + // examines at a 3x3 section of the water array around the passed location + // and only drains water away from center location + int rowloc, colloc; + double ht_diff, center_water_elev, cell_water_elev, flow; + + for (rowloc=centerrow-1; rowloc<=centerrow+1; rowloc++){ + for (colloc=centercol-1; colloc<=centercol+1; colloc++){ + // make sure centre element is not included + if (((rowloc != centerrow) || (colloc != centercol)) && + (bigdem[rowloc][colloc] > missingvalue)){ + cell_water_elev = bigdem[rowloc][colloc] + bigwater[rowloc][colloc]; + center_water_elev = bigdem[centerrow][centercol] + bigwater[centerrow][centercol]; + ht_diff = center_water_elev - cell_water_elev; + if (ht_diff > 0) { + if (bigdem[centerrow][centercol] > cell_water_elev) { + flow = bigwater[centerrow][centercol]/8.0; + } + else + { + flow = ((bigdem[centerrow][centercol] - bigdem[rowloc][colloc]) + + (bigwater[centerrow][centercol] - bigwater[rowloc][colloc]))/8.0; + } + flow = min(max(flow, 0.0), bigwater[centerrow][centercol]); + bigwater[centerrow][centercol] = max(bigwater[centerrow][centercol] - flow, 0.0); + bigwater[rowloc][colloc] = bigwater[rowloc][colloc] + flow; + } + } + } + } +} + + +void runoffs(int centerrow, int centercol, double missingvalue) { + // examines at a 3x3 section of the water array around the passed location + // and only drains water away from center location + int rowloc, colloc; + double ht_diff, center_water_elev, cell_water_elev, flow; + + for (rowloc=centerrow-1; rowloc<=centerrow+1; rowloc++){ + for (colloc=centercol-1; colloc<=centercol+1; colloc++){ + // make sure centre element is not included + if (((rowloc != centerrow) || (colloc != centercol)) && + (bigdem[rowloc][colloc] > missingvalue)){ + ht_diff = (bigdem[centerrow][centercol] + bigwater[centerrow][centercol])- + (bigdem[rowloc][colloc] + bigwater[rowloc][colloc]); + if (ht_diff > 0) { + if (bigdem[centerrow][centercol] > (bigdem[rowloc][colloc]+bigwater[rowloc][colloc])) { + flow = bigwater[centerrow][centercol]/8.0; + } + else + { + flow = ((bigdem[centerrow][centercol] - bigdem[rowloc][colloc]) + + (bigwater[centerrow][centercol] - bigwater[rowloc][colloc]))/8.0; + flow = ht_diff/8.0; + } + flow = min(flow, bigwater[centerrow][centercol]); + bigwater[centerrow][centercol] = bigwater[centerrow][centercol] - flow; + bigwater[rowloc][colloc] = bigwater[rowloc][colloc] + flow; + } + } + } + } +} + + +void runoffd(int centerrow, int centercol, int drainrow, int draincol, double missingvalue) { + // examines at a 3x3 section of the water array around the passed location + // and only drains water away from center location + int rowloc, colloc; + double ht_diff, center_water_elev, cell_water_elev, flow; + for (rowloc=centerrow-1; rowloc<=centerrow+1; rowloc++){ + for (colloc=centercol-1; colloc<=centercol+1; colloc++){ + // make sure centre element is not included + if (((rowloc != centerrow) || (colloc != centercol)) && + (bigdem[rowloc][colloc] > missingvalue)){ + center_water_elev = bigdem[centerrow][centercol] + bigwater[centerrow][centercol]; + cell_water_elev = bigdem[rowloc][colloc] + bigwater[rowloc][colloc]; + // check for drain + if((colloc == draincol) && (rowloc == drainrow)){ + // drain all water in center cell and edge cell + totaldrain = totaldrain + bigwater[drainrow][draincol] + bigwater[centerrow][centercol]; + bigwater[drainrow][draincol] = 0.0; + bigwater[centerrow][centercol] = 0.0; + } + else + { + ht_diff = center_water_elev - cell_water_elev; + if (ht_diff > 0) { + if (bigdem[centerrow][centercol] > cell_water_elev) { + flow = bigwater[centerrow][centercol]/8.0; + } + else + { + flow = ((bigdem[centerrow][centercol] - bigdem[rowloc][colloc]) + + (bigwater[centerrow][centercol] - bigwater[rowloc][colloc]))/8.0; + } + flow = min(max(flow, 0.0), bigwater[centerrow][centercol]); + bigwater[centerrow][centercol] = max(bigwater[centerrow][centercol] - flow, 0.0); + bigwater[rowloc][colloc] = bigwater[rowloc][colloc] + flow; + } + } + } + } + } +}