[GRASS-SVN] r40536 - grass/trunk/lib/gpde
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jan 18 14:37:19 EST 2010
Author: huhabla
Date: 2010-01-18 14:37:18 -0500 (Mon, 18 Jan 2010)
New Revision: 40536
Modified:
grass/trunk/lib/gpde/N_gwflow.c
grass/trunk/lib/gpde/N_gwflow.h
Log:
Water budget for 3d groundwater flow.
River and drainage bugfix.
Modified: grass/trunk/lib/gpde/N_gwflow.c
===================================================================
--- grass/trunk/lib/gpde/N_gwflow.c 2010-01-18 18:00:25 UTC (rev 40535)
+++ grass/trunk/lib/gpde/N_gwflow.c 2010-01-18 19:37:18 UTC (rev 40536)
@@ -348,7 +348,101 @@
return mat_pos;
}
+
/* *************************************************************** */
+/* ****************** N_gwflow_3d_calc_water_budget ************** */
+/* *************************************************************** */
+/*!
+ * \brief This function computes the water budget of the entire groundwater
+ *
+ * The water budget is calculated for each active and dirichlet cell from
+ * its surrounding neighbours. This is based on the 7 star mass balance computation
+ * of N_callback_gwflow_3d and the gradient of the water heights in the cells.
+ * The sum of the water budget of each active/dirichlet cell must be near zero
+ * due the effect of numerical inaccuracy of cpu's.
+ *
+ * \param gwdata N_gwflow_data3d *
+ * \param geom N_geom_data *
+ * \param budget N_array_3d
+ * \returnvoid
+ *
+ * */
+void
+N_gwflow_3d_calc_water_budget(N_gwflow_data3d * data, N_geom_data * geom, N_array_3d * budget)
+{
+ int z, y, x, stat;
+ double h, hc;
+ double val;
+ double sum;
+ N_data_star *dstar;
+
+ int rows = data->status->rows;
+ int cols = data->status->cols;
+ int depths = data->status->depths;
+ sum = 0;
+
+ for (z = 0; z < depths; z++) {
+ for (y = 0; y < rows; y++) {
+ G_percent(y, rows - 1, 10);
+ for (x = 0; x < cols; x++) {
+ stat = (int)N_get_array_3d_d_value(data->status, x, y);
+
+ val = 0.0;
+
+ if (stat != N_CELL_INACTIVE ) { /*all active/dirichlet cells */
+
+ /* Compute the flow parameter */
+ dstar = N_callback_gwflow_3d(data, geom, x, y, z);
+ /* Compute the gradient in each direction pointing from the center */
+ hc = N_get_array_3d_d_value(data->phead, x, y);
+
+ if((int)N_get_array_3d_d_value(data->status, x + 1, y, z ) != N_CELL_INACTIVE) {
+ h = N_get_array_3d_d_value(data->phead, x + 1, y, z );
+ val += dstar->E * (hc - h);
+ }
+ if((int)N_get_array_3d_d_value(data->status, x - 1, y, z ) != N_CELL_INACTIVE) {
+ h = N_get_array_3d_d_value(data->phead, x - 1, y, z);
+ val += dstar->W * (hc - h);
+ }
+ if((int)N_get_array_3d_d_value(data->status, x , y + 1, z) != N_CELL_INACTIVE) {
+ h = N_get_array_3d_d_value(data->phead, x , y + 1, z);
+ val += dstar->S * (hc - h);
+ }
+ if((int)N_get_array_3d_d_value(data->status, x , y - 1) != N_CELL_INACTIVE) {
+ h = N_get_array_3d_d_value(data->phead, x , y - 1);
+ val += dstar->N * (hc - h);
+ }
+ if((int)N_get_array_3d_d_value(data->status, x , y, z + 1) != N_CELL_INACTIVE) {
+ h = N_get_array_3d_d_value(data->phead, x , y, z + 1);
+ val += dstar->T * (hc - h);
+ }
+ if((int)N_get_array_3d_d_value(data->status, x , y, z - 1) != N_CELL_INACTIVE) {
+ h = N_get_array_3d_d_value(data->phead, x , y, z - 1);
+ val += dstar->B * (hc - h);
+ }
+ sum += val;
+
+ G_free(dstar);
+ }
+ else {
+ Rast_set_null_value(&val, 1, DCELL_TYPE);
+ }
+ N_put_array_3d_d_value(budget, x, y, z, val);
+ }
+ }
+ }
+
+ if(fabs(sum) < 0.0000000001)
+ G_message("The total sum of the water budget: %g\n", sum);
+ else
+ G_warning("The total sum of the water budget is significant larger then 0: %g\n", sum);
+
+ return;
+}
+
+
+
+/* *************************************************************** */
/* ****************** N_callback_gwflow_2d *********************** */
/* *************************************************************** */
/*!
@@ -480,9 +574,12 @@
T_n = N_calc_harmonic_mean(hc_yn, hc_y) * z_n;
T_s = N_calc_harmonic_mean(hc_ys, hc_y) * z_s;
- /*compute the river leakage, this is an explicit method */
+ /* Compute the river leakage, this is an explicit method
+ * Rivers are only enabled, if the river bed is lower or equal to the surface
+ */
if (data->river_leak &&
- (N_get_array_2d_d_value(data->river_leak, col, row) != 0)) {
+ (N_get_array_2d_d_value(data->river_leak, col, row) != 0) &&
+ N_get_array_2d_d_value(data->river_bed, col, row) <= top) {
if (hc > N_get_array_2d_d_value(data->river_bed, col, row)) {
river_vect = N_get_array_2d_d_value(data->river_head, col, row) *
N_get_array_2d_d_value(data->river_leak, col, row);
@@ -496,9 +593,12 @@
}
}
- /*compute the drainage, this is an explicit method */
+ /* compute the drainage, this is an explicit method
+ * Drainage is only enabled, if the drain bed is lower or equal to the surface
+ */
if (data->drain_leak &&
- (N_get_array_2d_d_value(data->drain_leak, col, row) != 0)) {
+ (N_get_array_2d_d_value(data->drain_leak, col, row) != 0) &&
+ N_get_array_2d_d_value(data->drain_bed, col, row) <= top) {
if (hc > N_get_array_2d_d_value(data->drain_bed, col, row)) {
drain_vect = N_get_array_2d_d_value(data->drain_bed, col, row) *
N_get_array_2d_d_value(data->drain_leak, col, row);
Modified: grass/trunk/lib/gpde/N_gwflow.h
===================================================================
--- grass/trunk/lib/gpde/N_gwflow.h 2010-01-18 18:00:25 UTC (rev 40535)
+++ grass/trunk/lib/gpde/N_gwflow.h 2010-01-18 19:37:18 UTC (rev 40536)
@@ -100,6 +100,8 @@
int col, int row, int depth);
extern N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom,
int col, int row);
+extern void N_gwflow_3d_calc_water_budget(N_gwflow_data3d * data,
+ N_geom_data * geom, N_array_3d * budget);
extern void N_gwflow_2d_calc_water_budget(N_gwflow_data2d * data,
N_geom_data * geom, N_array_2d * balance);
extern N_gwflow_data3d *N_alloc_gwflow_data3d(int cols, int rows, int depths,
More information about the grass-commit
mailing list