Loading...
Searching...
No Matches
n_gwflow.c
Go to the documentation of this file.
/*****************************************************************************
*
* MODULE: Grass PDE Numerical Library
* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
* soerengebbert <at> gmx <dot> de
*
* PURPOSE: groundwater flow in porous media
* part of the gpde library
*
* COPYRIGHT: (C) 2000 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*
*****************************************************************************/
#include <grass/N_gwflow.h>
/* *************************************************************** */
/* ***************** N_gwflow_data3d ***************************** */
/* *************************************************************** */
/*!
* \brief Allocate memory for the groundwater calculation data structure in 3
* dimensions
*
* The groundwater calculation data structure will be allocated including
* all appendant 3d and 2d arrays. The offset for the 3d arrays is one
* to establish homogeneous Neumann boundary conditions at the calculation area
* border. This data structure is used to create a linear equation system based
* on the computation of groundwater flow in porous media with the finite volume
* method.
*
* \param cols int
* \param rows int
* \param depths int
* \return N_gwflow_data3d *
* */
N_gwflow_data3d *N_alloc_gwflow_data3d(int cols, int rows, int depths,
int river, int drain)
{
N_gwflow_data3d *data;
data = (N_gwflow_data3d *)G_calloc(1, sizeof(N_gwflow_data3d));
data->phead = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->phead_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->hc_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->hc_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->hc_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->s = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
if (river) {
data->river_head = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->river_leak = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->river_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
}
else {
data->river_head = NULL;
data->river_leak = NULL;
data->river_bed = NULL;
}
if (drain) {
data->drain_leak = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->drain_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
}
else {
data->drain_leak = NULL;
data->drain_bed = NULL;
}
return data;
}
/* *************************************************************** */
/* ********************* N_free_gwflow_data3d ******************** */
/* *************************************************************** */
/*!
* \brief Release the memory of the groundwater flow data structure in three
* dimensions
*
* \param data N_gwflow_data3d *
* \return void *
* */
void N_free_gwflow_data3d(N_gwflow_data3d *data)
{
if (data->phead)
N_free_array_3d(data->phead);
if (data->phead_start)
N_free_array_3d(data->phead_start);
if (data->status)
N_free_array_3d(data->status);
if (data->hc_x)
N_free_array_3d(data->hc_x);
if (data->hc_y)
N_free_array_3d(data->hc_y);
if (data->hc_z)
N_free_array_3d(data->hc_z);
if (data->q)
N_free_array_3d(data->q);
if (data->s)
N_free_array_3d(data->s);
if (data->nf)
N_free_array_3d(data->nf);
if (data->r)
N_free_array_2d(data->r);
if (data->river_head)
N_free_array_3d(data->river_head);
if (data->river_leak)
N_free_array_3d(data->river_leak);
if (data->river_bed)
N_free_array_3d(data->river_bed);
if (data->drain_leak)
N_free_array_3d(data->drain_leak);
if (data->drain_bed)
N_free_array_3d(data->drain_bed);
G_free(data);
data = NULL;
return;
}
/* *************************************************************** */
/* ******************** N_alloc_gwflow_data2d ******************** */
/* *************************************************************** */
/*!
* \brief Allocate memory for the groundwater calculation data structure in 2
* dimensions
*
* The groundwater calculation data structure will be allocated including
* all appendant 2d arrays. The offset for the 3d arrays is one
* to establish homogeneous Neumann boundary conditions at the calculation area
* border. This data structure is used to create a linear equation system based
* on the computation of groundwater flow in porous media with the finite volume
* method.
*
* \param cols int
* \param rows int
* \param river
* \param drain
* \return N_gwflow_data2d *
* */
N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river, int drain)
{
N_gwflow_data2d *data;
data = (N_gwflow_data2d *)G_calloc(1, sizeof(N_gwflow_data2d));
data->phead = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->phead_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->status = N_alloc_array_2d(cols, rows, 1, CELL_TYPE);
data->hc_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->hc_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->s = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
if (river) {
data->river_head = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->river_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->river_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
}
else {
data->river_head = NULL;
data->river_leak = NULL;
data->river_bed = NULL;
}
if (drain) {
data->drain_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->drain_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
}
else {
data->drain_leak = NULL;
data->drain_bed = NULL;
}
return data;
}
/* *************************************************************** */
/* ****************** N_free_gwflow_data2d *********************** */
/* *************************************************************** */
/*!
* \brief Release the memory of the groundwater flow data structure in two
* dimensions
*
* \param data N_gwflow_data2d *
* \return void
* */
void N_free_gwflow_data2d(N_gwflow_data2d *data)
{
if (data->phead)
N_free_array_2d(data->phead);
if (data->phead_start)
N_free_array_2d(data->phead_start);
if (data->status)
N_free_array_2d(data->status);
if (data->hc_x)
N_free_array_2d(data->hc_x);
if (data->hc_y)
N_free_array_2d(data->hc_y);
if (data->q)
N_free_array_2d(data->q);
if (data->s)
N_free_array_2d(data->s);
if (data->nf)
N_free_array_2d(data->nf);
if (data->r)
N_free_array_2d(data->r);
if (data->top)
N_free_array_2d(data->top);
if (data->bottom)
N_free_array_2d(data->bottom);
if (data->river_head)
N_free_array_2d(data->river_head);
if (data->river_leak)
N_free_array_2d(data->river_leak);
if (data->river_bed)
N_free_array_2d(data->river_bed);
if (data->drain_leak)
N_free_array_2d(data->drain_leak);
if (data->drain_bed)
N_free_array_2d(data->drain_bed);
G_free(data);
data = NULL;
;
return;
}
/* *************************************************************** */
/* ***************** N_callback_gwflow_3d ************************ */
/* *************************************************************** */
/*!
* \brief This callback function creates the mass balance of a 7 point star
*
* The mass balance is based on the common groundwater flow equation:
*
* \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
*
* This equation is discretizised with the finite volume method in three
* dimensions.
*
*
* \param gwdata N_gwflow_data3d *
* \param geom N_geom_data *
* \param col int
* \param row int
* \param depth int
* \return N_data_star *
*
* */
N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data *geom, int col,
int row, int depth)
{
double hc_e = 0, hc_w = 0, hc_n = 0, hc_s = 0, hc_t = 0, hc_b = 0;
double dx, dy, dz, Ax, Ay, Az;
double hc_x, hc_y, hc_z;
double hc_xw, hc_yn, hc_zt;
double hc_xe, hc_ys, hc_zb;
double hc_start;
double Ss, r, /* nf, */ q;
double C, W, E, N, S, T, B, V;
N_data_star *mat_pos;
N_gwflow_data3d *data;
/*cast the void pointer to the right data structure */
data = (N_gwflow_data3d *)gwdata;
dx = geom->dx;
dy = geom->dy;
dz = geom->dz;
Az = N_get_geom_data_area_of_cell(geom, row);
Ay = geom->dx * geom->dz;
Ax = geom->dz * geom->dy;
/*read the data from the arrays */
hc_start = N_get_array_3d_d_value(data->phead_start, col, row, depth);
hc_x = N_get_array_3d_d_value(data->hc_x, col, row, depth);
hc_y = N_get_array_3d_d_value(data->hc_y, col, row, depth);
hc_z = N_get_array_3d_d_value(data->hc_z, col, row, depth);
hc_xw = N_get_array_3d_d_value(data->hc_x, col - 1, row, depth);
hc_xe = N_get_array_3d_d_value(data->hc_x, col + 1, row, depth);
hc_yn = N_get_array_3d_d_value(data->hc_y, col, row - 1, depth);
hc_ys = N_get_array_3d_d_value(data->hc_y, col, row + 1, depth);
hc_zt = N_get_array_3d_d_value(data->hc_z, col, row, depth + 1);
hc_zb = N_get_array_3d_d_value(data->hc_z, col, row, depth - 1);
hc_w = N_calc_harmonic_mean(hc_xw, hc_x);
hc_e = N_calc_harmonic_mean(hc_xe, hc_x);
hc_n = N_calc_harmonic_mean(hc_yn, hc_y);
hc_s = N_calc_harmonic_mean(hc_ys, hc_y);
hc_t = N_calc_harmonic_mean(hc_zt, hc_z);
hc_b = N_calc_harmonic_mean(hc_zb, hc_z);
/*inner sources */
q = N_get_array_3d_d_value(data->q, col, row, depth);
/*storativity */
Ss = N_get_array_3d_d_value(data->s, col, row, depth);
/*porosity */
/* nf = N_get_array_3d_d_value(data->nf, col, row, depth); */
/*mass balance center cell to western cell */
W = -1 * Ax * hc_w / dx;
/*mass balance center cell to eastern cell */
E = -1 * Ax * hc_e / dx;
/*mass balance center cell to northern cell */
N = -1 * Ay * hc_n / dy;
/*mass balance center cell to southern cell */
S = -1 * Ay * hc_s / dy;
/*mass balance center cell to top cell */
T = -1 * Az * hc_t / dz;
/*mass balance center cell to bottom cell */
B = -1 * Az * hc_b / dz;
/*storativity */
Ss = Az * dz * Ss;
/*the diagonal entry of the matrix */
C = -1 * (W + E + N + S + T + B - Ss / data->dt * Az);
/*the entry in the right side b of Ax = b */
V = (q + hc_start * Ss / data->dt * Az);
/*only the top cells will have recharge */
if (depth == geom->depths - 2) {
r = N_get_array_2d_d_value(data->r, col, row);
V += r * Az;
}
G_debug(5, "N_callback_gwflow_3d: called [%i][%i][%i]", depth, col, row);
/*create the 7 point star entries */
mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
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
* \return void
*
* */
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, z);
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, z);
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,
z) != N_CELL_INACTIVE) {
h = N_get_array_3d_d_value(data->phead, x, y - 1, z);
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 significantly larger "
"then 0: %g\n"),
sum);
return;
}
/* *************************************************************** */
/* ****************** N_callback_gwflow_2d *********************** */
/* *************************************************************** */
/*!
* \brief This callback function creates the mass balance of a 5 point star
*
* The mass balance is based on the common groundwater flow equation:
*
* \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
*
* This equation is discretizised with the finite volume method in two
* dimensions.
*
* \param gwdata N_gwflow_data2d *
* \param geom N_geom_data *
* \param col int
* \param row int
* \return N_data_star *
*
* */
N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data *geom, int col,
int row)
{
double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
double dx, dy, Az;
double hc_x, hc_y;
double z, top;
double hc_xw, hc_yn;
double z_xw, z_yn;
double hc_xe, hc_ys;
double z_xe, z_ys;
double hc, hc_start;
double Ss, r, q;
double C, W, E, N, S, V;
N_gwflow_data2d *data;
N_data_star *mat_pos;
double river_vect = 0; /*entry in vector */
double river_mat = 0; /*entry in matrix */
double drain_vect = 0; /*entry in vector */
double drain_mat = 0; /*entry in matrix */
/*cast the void pointer to the right data structure */
data = (N_gwflow_data2d *)gwdata;
dx = geom->dx;
dy = geom->dy;
Az = N_get_geom_data_area_of_cell(geom, row);
/*read the data from the arrays */
hc_start = N_get_array_2d_d_value(data->phead_start, col, row);
hc = N_get_array_2d_d_value(data->phead, col, row);
top = N_get_array_2d_d_value(data->top, col, row);
/* Inner sources */
q = N_get_array_2d_d_value(data->q, col, row);
/* storativity or porosity of current cell face [-] */
Ss = N_get_array_2d_d_value(data->s, col, row);
/* recharge */
r = N_get_array_2d_d_value(data->r, col, row) * Az;
if (hc > top) { /*If the aquifer is confined */
z = N_get_array_2d_d_value(data->top, col, row) -
N_get_array_2d_d_value(data->bottom, col, row);
z_xw = N_get_array_2d_d_value(data->top, col - 1, row) -
N_get_array_2d_d_value(data->bottom, col - 1, row);
z_xe = N_get_array_2d_d_value(data->top, col + 1, row) -
N_get_array_2d_d_value(data->bottom, col + 1, row);
z_yn = N_get_array_2d_d_value(data->top, col, row - 1) -
N_get_array_2d_d_value(data->bottom, col, row - 1);
z_ys = N_get_array_2d_d_value(data->top, col, row + 1) -
N_get_array_2d_d_value(data->bottom, col, row + 1);
}
else { /* the aquifer is unconfined */
/* If the aquifer is unconfied use an explicit scheme to solve
* the nonlinear equation. We use the phead from the first iteration */
z = N_get_array_2d_d_value(data->phead, col, row) -
N_get_array_2d_d_value(data->bottom, col, row);
z_xw = N_get_array_2d_d_value(data->phead, col - 1, row) -
N_get_array_2d_d_value(data->bottom, col - 1, row);
z_xe = N_get_array_2d_d_value(data->phead, col + 1, row) -
N_get_array_2d_d_value(data->bottom, col + 1, row);
z_yn = N_get_array_2d_d_value(data->phead, col, row - 1) -
N_get_array_2d_d_value(data->bottom, col, row - 1);
z_ys = N_get_array_2d_d_value(data->phead, col, row + 1) -
N_get_array_2d_d_value(data->bottom, col, row + 1);
}
/*geometrical mean of cell height */
if (z_w > 0 || z_w < 0 || z_w == 0)
z_w = N_calc_arith_mean(z_xw, z);
else
z_w = z;
if (z_e > 0 || z_e < 0 || z_e == 0)
z_e = N_calc_arith_mean(z_xe, z);
else
z_e = z;
if (z_n > 0 || z_n < 0 || z_n == 0)
z_n = N_calc_arith_mean(z_yn, z);
else
z_n = z;
if (z_s > 0 || z_s < 0 || z_s == 0)
z_s = N_calc_arith_mean(z_ys, z);
else
z_s = z;
/*get the surrounding permeabilities */
hc_x = N_get_array_2d_d_value(data->hc_x, col, row);
hc_y = N_get_array_2d_d_value(data->hc_y, col, row);
hc_xw = N_get_array_2d_d_value(data->hc_x, col - 1, row);
hc_xe = N_get_array_2d_d_value(data->hc_x, col + 1, row);
hc_yn = N_get_array_2d_d_value(data->hc_y, col, row - 1);
hc_ys = N_get_array_2d_d_value(data->hc_y, col, row + 1);
/* calculate the transmissivities */
T_w = N_calc_harmonic_mean(hc_xw, hc_x) * z_w;
T_e = N_calc_harmonic_mean(hc_xe, hc_x) * z_e;
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
* Influent and effluent flow is computed.
*/
if (data->river_leak &&
(N_get_array_2d_d_value(data->river_leak, col, row) != 0) &&
N_get_array_2d_d_value(data->river_bed, col, row) <= top) {
/* Groundwater surface is above the river bed */
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);
river_mat = N_get_array_2d_d_value(data->river_leak, col, row);
} /* Groundwater surface is below the river bed */
else 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_bed, col, row)) *
N_get_array_2d_d_value(data->river_leak, col, row);
river_mat = 0;
}
}
/* compute the drainage, this is an explicit method
* Drainage is only enabled, if the drain bed is lower the groundwater
* 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_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);
drain_mat = N_get_array_2d_d_value(data->drain_leak, col, row);
}
else if (hc <= N_get_array_2d_d_value(data->drain_bed, col, row)) {
drain_vect = 0;
drain_mat = 0;
}
}
/*mass balance center cell to western cell */
W = -1 * T_w * dy / dx;
/*mass balance center cell to eastern cell */
E = -1 * T_e * dy / dx;
/*mass balance center cell to northern cell */
N = -1 * T_n * dx / dy;
/*mass balance center cell to southern cell */
S = -1 * T_s * dx / dy;
/*the diagonal entry of the matrix */
C = -1 *
(W + E + N + S - Az * Ss / data->dt - river_mat * Az - drain_mat * Az);
/*the entry in the right side b of Ax = b */
V = (q + hc_start * Az * Ss / data->dt) + r + river_vect * Az +
drain_vect * Az;
G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
/*create the 5 point star entries */
mat_pos = N_create_5star(C, W, E, N, S, V);
return mat_pos;
}
/* *************************************************************** */
/* ****************** N_gwflow_2d_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 5 star mass balance
* computation of N_callback_gwflow_2d 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 data N_gwflow_data2d *
* \param geom N_geom_data *
* \param budget N_array_2d
* \return void
*
* */
void N_gwflow_2d_calc_water_budget(N_gwflow_data2d *data, N_geom_data *geom,
N_array_2d *budget)
{
int y, x, stat;
double h, hc;
double val;
double sum;
N_data_star *dstar;
int rows = data->status->rows;
int cols = data->status->cols;
sum = 0;
for (y = 0; y < rows; y++) {
G_percent(y, rows - 1, 10);
for (x = 0; x < cols; x++) {
stat = N_get_array_2d_c_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_2d(data, geom, x, y);
/* Compute the gradient in each direction pointing from the
* center */
hc = N_get_array_2d_d_value(data->phead, x, y);
if ((int)N_get_array_2d_d_value(data->status, x + 1, y) !=
N_CELL_INACTIVE) {
h = N_get_array_2d_d_value(data->phead, x + 1, y);
val += dstar->E * (hc - h);
}
if ((int)N_get_array_2d_d_value(data->status, x - 1, y) !=
N_CELL_INACTIVE) {
h = N_get_array_2d_d_value(data->phead, x - 1, y);
val += dstar->W * (hc - h);
}
if ((int)N_get_array_2d_d_value(data->status, x, y + 1) !=
N_CELL_INACTIVE) {
h = N_get_array_2d_d_value(data->phead, x, y + 1);
val += dstar->S * (hc - h);
}
if ((int)N_get_array_2d_d_value(data->status, x, y - 1) !=
N_CELL_INACTIVE) {
h = N_get_array_2d_d_value(data->phead, x, y - 1);
val += dstar->N * (hc - h);
}
sum += val;
G_free(dstar);
}
else {
Rast_set_null_value(&val, 1, DCELL_TYPE);
}
N_put_array_2d_d_value(budget, x, y, 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 significantly larger "
"then 0: %g\n"),
sum);
return;
}
N_CELL_INACTIVE
#define N_CELL_INACTIVE
Definition N_pde.h:30
N_calc_arith_mean
double N_calc_arith_mean(double a, double b)
Calculate the arithmetic mean of values a and b.
Definition n_tools.c:31
N_calc_harmonic_mean
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition n_tools.c:115
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
NULL
#define NULL
Definition ccmath.h:32
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
r
double r
Definition driver/set_window.c:5
G_message
void G_message(const char *msg,...)
Print a message to stderr.
Definition gis/error.c:89
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
N_alloc_array_3d
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure.
Definition n_arrays.c:719
N_get_array_2d_c_value
CELL N_get_array_2d_c_value(N_array_2d *data, int col, int row)
Returns the value of type CELL at position col, row.
Definition n_arrays.c:314
N_free_array_3d
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Definition n_arrays.c:774
N_get_array_2d_d_value
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
Definition n_arrays.c:380
N_free_array_2d
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
Definition n_arrays.c:132
N_put_array_3d_d_value
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth.
Definition n_arrays.c:1148
N_alloc_array_2d
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
Definition n_arrays.c:75
N_get_array_3d_d_value
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition n_arrays.c:979
N_put_array_2d_d_value
void N_put_array_2d_d_value(N_array_2d *data, int col, int row, DCELL value)
Writes a DCELL value to the N_array_2d struct at position col, row.
Definition n_arrays.c:576
N_get_geom_data_area_of_cell
double N_get_geom_data_area_of_cell(N_geom_data *geom, int row)
Get the areay size in square meter of one cell (x*y) at row.
Definition n_geom.c:196
N_callback_gwflow_2d
N_data_star * N_callback_gwflow_2d(void *gwdata, N_geom_data *geom, int col, int row)
This callback function creates the mass balance of a 5 point star.
Definition n_gwflow.c:476
N_gwflow_3d_calc_water_budget
void N_gwflow_3d_calc_water_budget(N_gwflow_data3d *data, N_geom_data *geom, N_array_3d *budget)
This function computes the water budget of the entire groundwater.
Definition n_gwflow.c:373
N_alloc_gwflow_data2d
N_gwflow_data2d * N_alloc_gwflow_data2d(int cols, int rows, int river, int drain)
Allocate memory for the groundwater calculation data structure in 2 dimensions.
Definition n_gwflow.c:151
N_gwflow_2d_calc_water_budget
void N_gwflow_2d_calc_water_budget(N_gwflow_data2d *data, N_geom_data *geom, N_array_2d *budget)
This function computes the water budget of the entire groundwater.
Definition n_gwflow.c:659
N_free_gwflow_data3d
void N_free_gwflow_data3d(N_gwflow_data3d *data)
Release the memory of the groundwater flow data structure in three dimensions.
Definition n_gwflow.c:91
N_alloc_gwflow_data3d
N_gwflow_data3d * N_alloc_gwflow_data3d(int cols, int rows, int depths, int river, int drain)
Allocate memory for the groundwater calculation data structure in 3 dimensions.
Definition n_gwflow.c:39
N_free_gwflow_data2d
void N_free_gwflow_data2d(N_gwflow_data2d *data)
Release the memory of the groundwater flow data structure in two dimensions.
Definition n_gwflow.c:202
N_callback_gwflow_3d
N_data_star * N_callback_gwflow_3d(void *gwdata, N_geom_data *geom, int col, int row, int depth)
This callback function creates the mass balance of a 7 point star.
Definition n_gwflow.c:267
N_create_5star
N_data_star * N_create_5star(double C, double W, double E, double N, double S, double V)
allocate and initialize a 5 point star data structure
Definition n_les_assemble.c:124
N_create_7star
N_data_star * N_create_7star(double C, double W, double E, double N, double S, double T, double B, double V)
allocate and initialize a 7 point star data structure
Definition n_les_assemble.c:159
G_percent
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition percent.c:61
N_array_2d
Definition N_pde.h:132
N_array_3d
Definition N_pde.h:175
N_data_star
Matrix entries for a mass balance 5/7/9 star system.
Definition N_pde.h:295
N_data_star::E
double E
Definition N_pde.h:298
N_data_star::S
double S
Definition N_pde.h:298
N_data_star::W
double W
Definition N_pde.h:298
N_data_star::T
double T
Definition N_pde.h:300
N_data_star::N
double N
Definition N_pde.h:298
N_data_star::B
double B
Definition N_pde.h:302
N_geom_data
Geometric information about the structured grid.
Definition N_pde.h:101
N_geom_data::dx
double dx
Definition N_pde.h:108
N_geom_data::depths
int depths
Definition N_pde.h:114
N_geom_data::dy
double dy
Definition N_pde.h:109
N_geom_data::dz
double dz
Definition N_pde.h:110
N_gwflow_data2d
This data structure contains all data needed to compute the groundwater mass balance in two dimension...
Definition N_gwflow.h:65
N_gwflow_data2d::phead
N_array_2d * phead
Definition N_gwflow.h:66
N_gwflow_data3d
This data structure contains all data needed to compute the groundwater mass balance in three dimensi...
Definition N_gwflow.h:34
N_gwflow_data3d::phead
N_array_3d * phead
Definition N_gwflow.h:35
x
#define x