Loading...
Searching...
No Matches
n_solute_transport.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: solute transport in porous media
* part of the gpde library
*
* COPYRIGHT: (C) 2007 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 <math.h>
#include <grass/N_solute_transport.h>
/* ************************************************************************* *
* ************************************************************************* *
* ************************************************************************* */
/*! \brief This is just a placeholder
*
* */
N_data_star *N_callback_solute_transport_3d(void *solutedata, N_geom_data *geom,
int col, int row, int depth)
{
double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0, Df_t = 0, Df_b = 0;
double dx, dy, dz, Az;
double diff_x, diff_y, diff_z;
double diff_xw, diff_yn;
double diff_xe, diff_ys;
double diff_zt, diff_zb;
double cin = 0, /* cg, */ cg_start;
double R, nf, cs, q;
double C, W, E, N, S, T, B, V;
double vw = 0, ve = 0, vn = 0, vs = 0, vt = 0, vb = 0;
double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0, Ds_t = 0, Ds_b = 0;
double Dw = 0, De = 0, Dn = 0, Ds = 0, Dt = 0, Db = 0;
double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5, rt = 0.5, rb = 0.5;
N_solute_transport_data3d *data = NULL;
N_data_star *mat_pos;
N_gradient_3d grad;
/*cast the void pointer to the right data structure */
data = (N_solute_transport_data3d *)solutedata;
N_get_gradient_3d(data->grad, &grad, col, row, depth);
dx = geom->dx;
dy = geom->dy;
dz = geom->dz;
Az = N_get_geom_data_area_of_cell(geom, row);
/*read the data from the arrays */
cg_start = N_get_array_3d_d_value(data->c_start, col, row, depth);
/* cg = N_get_array_3d_d_value(data->c, col, row, depth); */
/*get the surrounding diffusion tensor entries */
diff_x = N_get_array_3d_d_value(data->diff_x, col, row, depth);
diff_y = N_get_array_3d_d_value(data->diff_y, col, row, depth);
diff_z = N_get_array_3d_d_value(data->diff_z, col, row, depth);
diff_xw = N_get_array_3d_d_value(data->diff_x, col - 1, row, depth);
diff_xe = N_get_array_3d_d_value(data->diff_x, col + 1, row, depth);
diff_yn = N_get_array_3d_d_value(data->diff_y, col, row - 1, depth);
diff_ys = N_get_array_3d_d_value(data->diff_y, col, row + 1, depth);
diff_zt = N_get_array_3d_d_value(data->diff_z, col, row, depth + 1);
diff_zb = N_get_array_3d_d_value(data->diff_z, col, row, depth - 1);
/* calculate the diffusion on the cell borders using the harmonical mean */
Df_w = N_calc_harmonic_mean(diff_xw, diff_x);
Df_e = N_calc_harmonic_mean(diff_xe, diff_x);
Df_n = N_calc_harmonic_mean(diff_yn, diff_y);
Df_s = N_calc_harmonic_mean(diff_ys, diff_y);
Df_t = N_calc_harmonic_mean(diff_zt, diff_z);
Df_b = N_calc_harmonic_mean(diff_zb, diff_z);
/* calculate the dispersion */
/*todo */
/* calculate the velocity parts with full upwinding scheme */
vw = grad.WC;
ve = grad.EC;
vn = grad.NC;
vs = grad.SC;
vt = grad.TC;
vb = grad.BC;
/* put the diffusion and dispersion together */
Dw = ((Df_w + Ds_w)) / dx;
De = ((Df_e + Ds_e)) / dx;
Dn = ((Df_n + Ds_n)) / dy;
Ds = ((Df_s + Ds_s)) / dy;
Dt = ((Df_t + Ds_t)) / dz;
Db = ((Df_b + Ds_b)) / dz;
rw = N_exp_upwinding(-1 * vw, dx, Dw);
re = N_exp_upwinding(ve, dx, De);
rs = N_exp_upwinding(-1 * vs, dy, Ds);
rn = N_exp_upwinding(vn, dy, Dn);
rb = N_exp_upwinding(-1 * vb, dz, Dn);
rt = N_exp_upwinding(vt, dz, Dn);
/*mass balance center cell to western cell */
W = -1 * (Dw)*dy * dz - vw * (1 - rw) * dy * dz;
/*mass balance center cell to eastern cell */
E = -1 * (De)*dy * dz + ve * (1 - re) * dy * dz;
/*mass balance center cell to southern cell */
S = -1 * (Ds)*dx * dz - vs * (1 - rs) * dx * dz;
/*mass balance center cell to northern cell */
N = -1 * (Dn)*dx * dz + vn * (1 - rn) * dx * dz;
/*mass balance center cell to bottom cell */
B = -1 * (Db)*Az - vb * (1 - rb) * Az;
/*mass balance center cell to top cell */
T = -1 * (Dt)*Az + vt * (1 - rt) * Az;
/* Retardation */
R = N_get_array_3d_d_value(data->R, col, row, depth);
/* Inner sources */
cs = N_get_array_3d_d_value(data->cs, col, row, depth);
/* effective porosity */
nf = N_get_array_3d_d_value(data->nf, col, row, depth);
/* groundwater sources and sinks */
q = N_get_array_3d_d_value(data->q, col, row, depth);
/* concentration of influent water */
cin = N_get_array_3d_d_value(data->cin, col, row, depth);
/*the diagonal entry of the matrix */
C = ((Dw - vw) * dy * dz + (De + ve) * dy * dz + (Ds - vs) * dx * dz +
(Dn + vn) * dx * dz + (Db - vb) * Az + (Dt + vt) * Az +
Az * dz * R / data->dt - q / nf);
/*the entry in the right side b of Ax = b */
V = (cs + cg_start * Az * dz * R / data->dt - q / nf * cin);
/*
* printf("nf %g\n", nf);
* printf("q %g\n", q);
* printf("cs %g\n", cs);
* printf("cin %g\n", cin);
* printf("cg %g\n", cg);
* printf("cg_start %g\n", cg_start);
* printf("Az %g\n", Az);
* printf("z %g\n", z);
* printf("R %g\n", R);
* printf("dt %g\n", data->dt);
*/
G_debug(6, "N_callback_solute_transport_3d: called [%i][%i][%i]", row, col,
depth);
/*create the 7 point star entries */
mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
return mat_pos;
}
/* ************************************************************************* *
* ************************************************************************* *
* ************************************************************************* */
/*!
* \brief This callback function creates the mass balance of a 5 point star
*
* The mass balance is based on the common solute transport equation:
*
* \f[\frac{\partial c_g}{\partial t} R = \nabla \cdot ({\bf D} \nabla c_g -
* {\bf u} c_g) + \sigma + \frac{q}{n_f}(c_g - c_in) \f]
*
* This equation is discretizised with the finite volume method in two
* dimensions.
*
*
* \param solutedata * N_solute_transport_data2d - a void pointer to the data
* structure \param geom N_geom_data * \param col int \param row int \return
* N_data_star * - a five point data star
*
* */
N_data_star *N_callback_solute_transport_2d(void *solutedata, N_geom_data *geom,
int col, int row)
{
double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0;
double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
double dx, dy, Az;
double diff_x, diff_y;
double disp_x, disp_y;
double z;
double diff_xw, diff_yn;
double disp_xw, disp_yn;
double z_xw, z_yn;
double diff_xe, diff_ys;
double disp_xe, disp_ys;
double z_xe, z_ys;
double cin = 0, /* cg, */ cg_start;
double R, nf, cs, q;
double C, W, E, N, S, V, NE, NW, SW, SE;
double vw = 0, ve = 0, vn = 0, vs = 0;
double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0;
double Dw = 0, De = 0, Dn = 0, Ds = 0;
double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5;
N_solute_transport_data2d *data = NULL;
N_data_star *mat_pos;
N_gradient_2d grad;
/*cast the void pointer to the right data structure */
data = (N_solute_transport_data2d *)solutedata;
N_get_gradient_2d(data->grad, &grad, col, row);
dx = geom->dx;
dy = geom->dy;
Az = N_get_geom_data_area_of_cell(geom, row);
/*read the data from the arrays */
cg_start = N_get_array_2d_d_value(data->c_start, col, row);
/* cg = N_get_array_2d_d_value(data->c, col, row); */
/* calculate the cell height */
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);
/*geometrical mean of cell height */
z_w = N_calc_geom_mean(z_xw, z);
z_e = N_calc_geom_mean(z_xe, z);
z_n = N_calc_geom_mean(z_yn, z);
z_s = N_calc_geom_mean(z_ys, z);
/*get the surrounding diffusion tensor entries */
diff_x = N_get_array_2d_d_value(data->diff_x, col, row);
diff_y = N_get_array_2d_d_value(data->diff_y, col, row);
diff_xw = N_get_array_2d_d_value(data->diff_x, col - 1, row);
diff_xe = N_get_array_2d_d_value(data->diff_x, col + 1, row);
diff_yn = N_get_array_2d_d_value(data->diff_y, col, row - 1);
diff_ys = N_get_array_2d_d_value(data->diff_y, col, row + 1);
/* calculate the diffusion at the cell borders using the harmonical mean */
Df_w = N_calc_harmonic_mean(diff_xw, diff_x);
Df_e = N_calc_harmonic_mean(diff_xe, diff_x);
Df_n = N_calc_harmonic_mean(diff_yn, diff_y);
Df_s = N_calc_harmonic_mean(diff_ys, diff_y);
/* calculate the dispersion */
/*get the surrounding dispersion tensor entries */
disp_x = N_get_array_2d_d_value(data->disp_xx, col, row);
disp_y = N_get_array_2d_d_value(data->disp_yy, col, row);
if (N_get_array_2d_d_value(data->status, col - 1, row) ==
N_CELL_TRANSMISSION) {
disp_xw = disp_x;
}
else {
disp_xw = N_get_array_2d_d_value(data->disp_xx, col - 1, row);
}
if (N_get_array_2d_d_value(data->status, col + 1, row) ==
N_CELL_TRANSMISSION) {
disp_xe = disp_x;
}
else {
disp_xe = N_get_array_2d_d_value(data->disp_xx, col + 1, row);
}
if (N_get_array_2d_d_value(data->status, col, row - 1) ==
N_CELL_TRANSMISSION) {
disp_yn = disp_y;
}
else {
disp_yn = N_get_array_2d_d_value(data->disp_yy, col, row - 1);
}
if (N_get_array_2d_d_value(data->status, col, row + 1) ==
N_CELL_TRANSMISSION) {
disp_ys = disp_y;
}
else {
disp_ys = N_get_array_2d_d_value(data->disp_yy, col, row + 1);
}
/* calculate the dispersion at the cell borders using the harmonical mean */
Ds_w = N_calc_harmonic_mean(disp_xw, disp_x);
Ds_e = N_calc_harmonic_mean(disp_xe, disp_x);
Ds_n = N_calc_harmonic_mean(disp_yn, disp_y);
Ds_s = N_calc_harmonic_mean(disp_ys, disp_y);
/* put the diffusion and dispersion together */
Dw = ((Df_w + Ds_w)) / dx;
De = ((Df_e + Ds_e)) / dx;
Ds = ((Df_s + Ds_s)) / dy;
Dn = ((Df_n + Ds_n)) / dy;
vw = -1.0 * grad.WC;
ve = grad.EC;
vs = -1.0 * grad.SC;
vn = grad.NC;
if (data->stab == N_UPWIND_FULL) {
rw = N_full_upwinding(vw, dx, Dw);
re = N_full_upwinding(ve, dx, De);
rs = N_full_upwinding(vs, dy, Ds);
rn = N_full_upwinding(vn, dy, Dn);
}
else if (data->stab == N_UPWIND_EXP) {
rw = N_exp_upwinding(vw, dx, Dw);
re = N_exp_upwinding(ve, dx, De);
rs = N_exp_upwinding(vs, dy, Ds);
rn = N_exp_upwinding(vn, dy, Dn);
}
/*mass balance center cell to western cell */
W = -1 * (Dw)*dy * z_w + vw * (1 - rw) * dy * z_w;
/*mass balance center cell to eastern cell */
E = -1 * (De)*dy * z_e + ve * (1 - re) * dy * z_e;
/*mass balance center cell to southern cell */
S = -1 * (Ds)*dx * z_s + vs * (1 - rs) * dx * z_s;
/*mass balance center cell to northern cell */
N = -1 * (Dn)*dx * z_n + vn * (1 - rn) * dx * z_n;
NW = 0.0;
SW = 0.0;
NE = 0.0;
SE = 0.0;
/* Retardation */
R = N_get_array_2d_d_value(data->R, col, row);
/* Inner sources */
cs = N_get_array_2d_d_value(data->cs, col, row);
/* effective porosity */
nf = N_get_array_2d_d_value(data->nf, col, row);
/* groundwater sources and sinks */
q = N_get_array_2d_d_value(data->q, col, row);
/* concentration of influent water */
cin = N_get_array_2d_d_value(data->cin, col, row);
/*the diagonal entry of the matrix */
C = (Dw + vw * rw) * dy * z_w + (De + ve * re) * dy * z_e +
(Ds + vs * rs) * dx * z_s + (Dn + vn * rn) * dx * z_n +
Az * z * R / data->dt - q / nf;
/*the entry in the right side b of Ax = b */
V = (cs + cg_start * Az * z * R / data->dt + q / nf * cin);
/*
fprintf(stderr, "nf %g\n", nf);
fprintf(stderr, "q %g\n", q);
fprintf(stderr, "cs %g\n", cs);
fprintf(stderr, "cin %g\n", cin);
fprintf(stderr, "cg %g\n", cg);
fprintf(stderr, "cg_start %g\n", cg_start);
fprintf(stderr, "Az %g\n", Az);
fprintf(stderr, "z %g\n", z);
fprintf(stderr, "R %g\n", R);
fprintf(stderr, "dt %g\n", data->dt);
*/
G_debug(6, "N_callback_solute_transport_2d: called [%i][%i]", row, col);
/*create the 9 point star entries */
mat_pos = N_create_9star(C, W, E, N, S, NW, SW, NE, SE, V);
return mat_pos;
}
/* ************************************************************************* *
* ************************************************************************* *
* ************************************************************************* */
/*!
* \brief Allocate memory for the solute transport data structure in three
* dimensions
*
* The solute transport data structure will be allocated including
* all appendant 3d 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 solute transport in porous media with the finite volume
* method.
*
* \param cols int
* \param rows int
* \param depths int
* \return N_solute_transport_data3d *
* */
N_solute_transport_data3d *N_alloc_solute_transport_data3d(int cols, int rows,
int depths)
{
N_solute_transport_data3d *data = NULL;
data = (N_solute_transport_data3d *)G_calloc(
1, sizeof(N_solute_transport_data3d));
data->c = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->c_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->diff_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->diff_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->diff_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->cs = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->R = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->cin = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
/*Allocate the dispersivity tensor */
data->disp_xx = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->disp_yy = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->disp_zz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->disp_xy = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->disp_xz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->disp_yz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
data->grad = N_alloc_gradient_field_3d(cols, rows, depths);
data->stab = N_UPWIND_EXP;
return data;
}
/* ************************************************************************* *
* ************************************************************************* *
* ************************************************************************* */
/*!
* \brief Allocate memory for the solute transport data structure in two
* dimensions
*
* The solute transport data structure will be allocated including
* all appendant 2d arrays. The offset for the 2d 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 solute transport in porous media with the finite volume
* method.
*
* \param cols int
* \param rows int
* \return N_solute_transport_data2d *
* */
N_solute_transport_data2d *N_alloc_solute_transport_data2d(int cols, int rows)
{
N_solute_transport_data2d *data = NULL;
data = (N_solute_transport_data2d *)G_calloc(
1, sizeof(N_solute_transport_data2d));
data->c = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->c_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->status = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->diff_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->diff_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->cs = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->R = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->cin = 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);
/*Allocate the dispersivity tensor */
data->disp_xx = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->disp_yy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->disp_xy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
data->grad = N_alloc_gradient_field_2d(cols, rows);
data->stab = N_UPWIND_EXP;
return data;
}
/* ************************************************************************* *
* ************************************************************************* *
* ************************************************************************* */
/*!
* \brief Release the memory of the solute transport data structure in three
* dimensions
*
* \param data N_solute_transport_data2d *
* \return void *
* */
void N_free_solute_transport_data3d(N_solute_transport_data3d *data)
{
N_free_array_3d(data->c);
N_free_array_3d(data->c_start);
N_free_array_3d(data->status);
N_free_array_3d(data->diff_x);
N_free_array_3d(data->diff_y);
N_free_array_3d(data->diff_z);
N_free_array_3d(data->q);
N_free_array_3d(data->cs);
N_free_array_3d(data->R);
N_free_array_3d(data->nf);
N_free_array_3d(data->cin);
N_free_array_3d(data->disp_xx);
N_free_array_3d(data->disp_yy);
N_free_array_3d(data->disp_zz);
N_free_array_3d(data->disp_xy);
N_free_array_3d(data->disp_xz);
N_free_array_3d(data->disp_yz);
G_free(data);
data = NULL;
return;
}
/* ************************************************************************* *
* ************************************************************************* *
* ************************************************************************* */
/*!
* \brief Release the memory of the solute transport data structure in two
* dimensions
*
* \param data N_solute_transport_data2d *
* \return void *
* */
void N_free_solute_transport_data2d(N_solute_transport_data2d *data)
{
N_free_array_2d(data->c);
N_free_array_2d(data->c_start);
N_free_array_2d(data->status);
N_free_array_2d(data->diff_x);
N_free_array_2d(data->diff_y);
N_free_array_2d(data->q);
N_free_array_2d(data->cs);
N_free_array_2d(data->R);
N_free_array_2d(data->nf);
N_free_array_2d(data->cin);
N_free_array_2d(data->top);
N_free_array_2d(data->bottom);
N_free_array_2d(data->disp_xx);
N_free_array_2d(data->disp_yy);
N_free_array_2d(data->disp_xy);
G_free(data);
data = NULL;
return;
}
/*!
* \brief Compute the transmission boundary condition in 2d
*
* This function calculates the transmission boundary condition
* for each cell with status N_CELL_TRANSMISSION. The surrounding
* gradient field is used to verify the flow direction. If a flow
* goes into a cell, the concentration (data->c) from the neighbour cell is
* added to the transmission cell. If the flow from several neighbour
* cells goes into the cell, the concentration mean is calculated.
*
* The new concentrations are written into the data->c_start array,
* so they can be handled by the matrix assembling function.
*
* \param data N_solute_transport_data2d *
* \return void *
* */
void N_calc_solute_transport_transmission_2d(N_solute_transport_data2d *data)
{
int i, j, count = 1;
int cols, rows;
double c;
N_gradient_2d grad;
cols = data->grad->cols;
rows = data->grad->rows;
G_debug(2, "N_calc_solute_transport_transmission_2d: calculating "
"transmission boundary");
for (j = 0; j < rows; j++) {
for (i = 0; i < cols; i++) {
if (N_get_array_2d_d_value(data->status, i, j) ==
N_CELL_TRANSMISSION) {
count = 0;
/*get the gradient neighbours */
N_get_gradient_2d(data->grad, &grad, i, j);
c = 0;
/*
c = N_get_array_2d_d_value(data->c_start, i, j);
if(c > 0)
count++;
*/
if (grad.WC > 0 &&
!N_is_array_2d_value_null(data->c, i - 1, j)) {
c += N_get_array_2d_d_value(data->c, i - 1, j);
count++;
}
if (grad.EC < 0 &&
!N_is_array_2d_value_null(data->c, i + 1, j)) {
c += N_get_array_2d_d_value(data->c, i + 1, j);
count++;
}
if (grad.NC < 0 &&
!N_is_array_2d_value_null(data->c, i, j - 1)) {
c += N_get_array_2d_d_value(data->c, i, j - 1);
count++;
}
if (grad.SC > 0 &&
!N_is_array_2d_value_null(data->c, i, j + 1)) {
c += N_get_array_2d_d_value(data->c, i, j + 1);
count++;
}
if (count != 0)
c = c / (double)count;
/*make sure it is not NAN */
if (c > 0 || c == 0 || c < 0)
N_put_array_2d_d_value(data->c_start, i, j, c);
}
}
}
return;
}
/*!
* \brief Compute the dispersivity tensor based on the solute transport data in
* 2d
*
* The dispersivity tensor is stored in the data structure.
* To compute the dispersivity tensor, the dispersivity lengths and the gradient
* field must be present.
*
* This is just a simple tensor computation which should be extended.
*
* \todo Change the tensor calculation to a more realistic algorithm
*
* \param data N_solute_transport_data2d *
* \return void *
* */
void N_calc_solute_transport_disptensor_2d(N_solute_transport_data2d *data)
{
int i, j;
int cols, rows;
double vx, vy, vv;
double disp_xx, disp_yy, disp_xy;
N_gradient_2d grad;
cols = data->grad->cols;
rows = data->grad->rows;
G_debug(2, "N_calc_solute_transport_disptensor_2d: calculating the "
"dispersivity tensor");
for (j = 0; j < rows; j++) {
for (i = 0; i < cols; i++) {
disp_xx = 0;
disp_yy = 0;
disp_xy = 0;
/*get the gradient neighbours */
N_get_gradient_2d(data->grad, &grad, i, j);
vx = (grad.WC + grad.EC) / 2;
vy = (grad.NC + grad.SC) / 2;
vv = sqrt(vx * vx + vy * vy);
if (vv != 0) {
disp_xx = data->al * vx * vx / vv + data->at * vy * vy / vv;
disp_yy = data->at * vx * vx / vv + data->al * vy * vy / vv;
disp_xy = (data->al - data->at) * vx * vy / vv;
}
G_debug(5,
"N_calc_solute_transport_disptensor_2d: [%i][%i] disp_xx "
"%g disp_yy %g disp_xy %g",
i, j, disp_xx, disp_yy, disp_xy);
N_put_array_2d_d_value(data->disp_xx, i, j, disp_xx);
N_put_array_2d_d_value(data->disp_yy, i, j, disp_yy);
N_put_array_2d_d_value(data->disp_xy, i, j, disp_xy);
}
}
return;
}
/*!
* \brief Compute the dispersivity tensor based on the solute transport data in
* 3d
*
* The dispersivity tensor is stored in the data structure.
* To compute the dispersivity tensor, the dispersivity lengths and the gradient
* field must be present.
*
* This is just a simple tensor computation which should be extended.
*
* \todo Change the tensor calculation to a more realistic algorithm
*
* \param data N_solute_transport_data3d *
* \return void *
* */
void N_calc_solute_transport_disptensor_3d(N_solute_transport_data3d *data)
{
int i, j, k;
int cols, rows, depths;
double vx, vy, vz, vv;
double disp_xx, disp_yy, disp_zz, disp_xy, disp_xz, disp_yz;
N_gradient_3d grad;
cols = data->grad->cols;
rows = data->grad->rows;
depths = data->grad->depths;
G_debug(2, "N_calc_solute_transport_disptensor_3d: calculating the "
"dispersivity tensor");
for (k = 0; k < depths; k++) {
for (j = 0; j < rows; j++) {
for (i = 0; i < cols; i++) {
disp_xx = 0;
disp_yy = 0;
disp_zz = 0;
disp_xy = 0;
disp_xz = 0;
disp_yz = 0;
/*get the gradient neighbours */
N_get_gradient_3d(data->grad, &grad, i, j, k);
vx = (grad.WC + grad.EC) / 2;
vy = (grad.NC + grad.SC) / 2;
vz = (grad.BC + grad.TC) / 2;
vv = sqrt(vx * vx + vy * vy + vz * vz);
if (vv != 0) {
disp_xx = data->al * vx * vx / vv +
data->at * vy * vy / vv + data->at * vz * vz / vv;
disp_yy = data->at * vx * vx / vv +
data->al * vy * vy / vv + data->at * vz * vz / vv;
disp_zz = data->at * vx * vx / vv +
data->at * vy * vy / vv + data->al * vz * vz / vv;
disp_xy = (data->al - data->at) * vx * vy / vv;
disp_xz = (data->al - data->at) * vx * vz / vv;
disp_yz = (data->al - data->at) * vy * vz / vv;
}
G_debug(5,
"N_calc_solute_transport_disptensor_3d: [%i][%i][%i] "
"disp_xx %g disp_yy %g disp_zz %g disp_xy %g disp_xz "
"%g disp_yz %g ",
i, j, k, disp_xx, disp_yy, disp_zz, disp_xy, disp_xz,
disp_yz);
N_put_array_3d_d_value(data->disp_xx, i, j, k, disp_xx);
N_put_array_3d_d_value(data->disp_yy, i, j, k, disp_yy);
N_put_array_3d_d_value(data->disp_zz, i, j, k, disp_zz);
N_put_array_3d_d_value(data->disp_xy, i, j, k, disp_xy);
N_put_array_3d_d_value(data->disp_xz, i, j, k, disp_xz);
N_put_array_3d_d_value(data->disp_yz, i, j, k, disp_yz);
}
}
}
return;
}
N_UPWIND_FULL
#define N_UPWIND_FULL
Definition N_pde.h:54
N_CELL_TRANSMISSION
#define N_CELL_TRANSMISSION
Definition N_pde.h:33
N_exp_upwinding
double N_exp_upwinding(double sprod, double distance, double D)
exponential upwinding stabilization algorithm
Definition n_upwind.c:63
N_calc_geom_mean
double N_calc_geom_mean(double a, double b)
Calculate the geometrical mean of values a and b.
Definition n_tools.c:73
N_full_upwinding
double N_full_upwinding(double sprod, double distance, double D)
full upwinding stabilization algorithm
Definition n_upwind.c:32
N_UPWIND_EXP
#define N_UPWIND_EXP
Definition N_pde.h:55
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
SE
#define SE
Definition dataquad.h:31
SW
#define SW
Definition dataquad.h:30
NE
#define NE
Definition dataquad.h:29
NW
#define NW
Definition dataquad.h:28
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
count
int count
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_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_is_array_2d_value_null
int N_is_array_2d_value_null(N_array_2d *data, int col, int row)
Returns 1 if the value of N_array_2d struct at position col, row is of type null, otherwise 0.
Definition n_arrays.c:231
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_alloc_gradient_field_2d
N_gradient_field_2d * N_alloc_gradient_field_2d(int cols, int rows)
Allocate a N_gradient_field_2d.
Definition n_gradient.c:896
N_alloc_gradient_field_3d
N_gradient_field_3d * N_alloc_gradient_field_3d(int cols, int rows, int depths)
Allocate a N_gradient_field_3d.
Definition n_gradient.c:993
N_get_gradient_2d
N_gradient_2d * N_get_gradient_2d(N_gradient_field_2d *field, N_gradient_2d *gradient, int col, int row)
Return a N_gradient_2d structure calculated from the input gradient field at position [row][col].
Definition n_gradient.c:114
N_get_gradient_3d
N_gradient_3d * N_get_gradient_3d(N_gradient_field_3d *field, N_gradient_3d *gradient, int col, int row, int depth)
Return a N_gradient_3d structure calculated from the input gradient field at position [depth][row][co...
Definition n_gradient.c:246
N_create_9star
N_data_star * N_create_9star(double C, double W, double E, double N, double S, double NW, double SW, double NE, double SE, double V)
allocate and initialize a 9 point star data structure
Definition n_les_assemble.c:200
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
N_alloc_solute_transport_data2d
N_solute_transport_data2d * N_alloc_solute_transport_data2d(int cols, int rows)
Allocate memory for the solute transport data structure in two dimensions.
Definition n_solute_transport.c:442
N_calc_solute_transport_disptensor_2d
void N_calc_solute_transport_disptensor_2d(N_solute_transport_data2d *data)
Compute the dispersivity tensor based on the solute transport data in 2d.
Definition n_solute_transport.c:637
N_free_solute_transport_data2d
void N_free_solute_transport_data2d(N_solute_transport_data2d *data)
Release the memory of the solute transport data structure in two dimensions.
Definition n_solute_transport.c:521
N_alloc_solute_transport_data3d
N_solute_transport_data3d * N_alloc_solute_transport_data3d(int cols, int rows, int depths)
Allocate memory for the solute transport data structure in three dimensions.
Definition n_solute_transport.c:389
N_free_solute_transport_data3d
void N_free_solute_transport_data3d(N_solute_transport_data3d *data)
Release the memory of the solute transport data structure in three dimensions.
Definition n_solute_transport.c:483
N_callback_solute_transport_2d
N_data_star * N_callback_solute_transport_2d(void *solutedata, N_geom_data *geom, int col, int row)
This callback function creates the mass balance of a 5 point star.
Definition n_solute_transport.c:180
N_callback_solute_transport_3d
N_data_star * N_callback_solute_transport_3d(void *solutedata, N_geom_data *geom, int col, int row, int depth)
This is just a placeholder.
Definition n_solute_transport.c:27
N_calc_solute_transport_transmission_2d
void N_calc_solute_transport_transmission_2d(N_solute_transport_data2d *data)
Compute the transmission boundary condition in 2d.
Definition n_solute_transport.c:563
N_calc_solute_transport_disptensor_3d
void N_calc_solute_transport_disptensor_3d(N_solute_transport_data3d *data)
Compute the dispersivity tensor based on the solute transport data in 3d.
Definition n_solute_transport.c:698
N_data_star
Matrix entries for a mass balance 5/7/9 star system.
Definition N_pde.h:295
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::dy
double dy
Definition N_pde.h:109
N_geom_data::dz
double dz
Definition N_pde.h:110
N_gradient_2d
Gradient between the cells in X and Y direction.
Definition N_pde.h:457
N_gradient_2d::EC
double EC
Definition N_pde.h:459
N_gradient_2d::SC
double SC
Definition N_pde.h:459
N_gradient_2d::WC
double WC
Definition N_pde.h:459
N_gradient_2d::NC
double NC
Definition N_pde.h:459
N_gradient_3d
Gradient between the cells in X, Y and Z direction.
Definition N_pde.h:464
N_gradient_3d::WC
double WC
Definition N_pde.h:466
N_gradient_3d::NC
double NC
Definition N_pde.h:466
N_gradient_3d::EC
double EC
Definition N_pde.h:466
N_gradient_3d::SC
double SC
Definition N_pde.h:466
N_gradient_3d::BC
double BC
Definition N_pde.h:466
N_gradient_3d::TC
double TC
Definition N_pde.h:466
N_solute_transport_data2d
Definition N_solute_transport.h:52
N_solute_transport_data3d
Definition N_solute_transport.h:23