19 #include "grass/N_solute_transport.h"
32 double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0, Df_t = 0, Df_b = 0;
34 double dx, dy, dz, Az;
36 double diff_x, diff_y, diff_z;
38 double diff_xw, diff_yn;
40 double diff_xe, diff_ys;
42 double diff_zt, diff_zb;
44 double cin = 0, cg, cg_start;
48 double C, W, E,
N, S, T, B, V;
50 double vw = 0, ve = 0, vn = 0, vs = 0, vt = 0, vb = 0;
52 double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0, Ds_t = 0, Ds_b = 0;
54 double Dw = 0, De = 0, Dn = 0, Ds = 0, Dt = 0, Db = 0;
56 double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5, rt = 0.5, rb = 0.5;
109 Dw = ((Df_w + Ds_w)) / dx;
110 De = ((Df_e + Ds_e)) / dx;
111 Dn = ((Df_n + Ds_n)) / dy;
112 Ds = ((Df_s + Ds_s)) / dy;
113 Dt = ((Df_t + Ds_t)) / dz;
114 Db = ((Df_b + Ds_b)) / dz;
124 W = -1 * (Dw) * dy * dz - vw * (1 - rw) * dy * dz;
126 E = -1 * (De) * dy * dz + ve * (1 - re) * dy * dz;
128 S = -1 * (Ds) * dx * dz - vs * (1 - rs) * dx * dz;
130 N = -1 * (Dn) * dx * dz + vn * (1 - rn) * dx * dz;
132 B = -1 * (Db) * Az - vb * (1 - rb) * Az;
134 T = -1 * (Dt) * Az + vt * (1 - rt) * Az;
148 C = ((Dw - vw) * dy * dz +
149 (De + ve) * dy * dz +
150 (Ds - vs) * dx * dz +
151 (Dn + vn) * dx * dz +
152 (Db - vb) * Az + (Dt + vt) * Az + Az * dz * R / data->
dt - q / nf);
155 V = (cs + cg_start * Az * dz * R / data->
dt - q / nf * cin);
169 G_debug(6,
"N_callback_solute_transport_3d: called [%i][%i][%i]", row,
202 double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0;
204 double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
208 double diff_x, diff_y;
210 double disp_x, disp_y;
214 double diff_xw, diff_yn;
216 double disp_xw, disp_yn;
220 double diff_xe, diff_ys;
222 double disp_xe, disp_ys;
226 double cin = 0, cg, cg_start;
232 double vw = 0, ve = 0, vn = 0, vs = 0;
234 double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0;
236 double Dw = 0, De = 0, Dn = 0, Ds = 0;
238 double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5;
340 Dw = ((Df_w + Ds_w)) / dx;
341 De = ((Df_e + Ds_e)) / dx;
342 Ds = ((Df_s + Ds_s)) / dy;
343 Dn = ((Df_n + Ds_n)) / dy;
364 W = -1 * (Dw) * dy * z_w + vw * (1 - rw) * dy * z_w;
366 E = -1 * (De) * dy * z_e + ve * (1 - re) * dy * z_e;
368 S = -1 * (Ds) * dx * z_s + vs * (1 - rs) * dx * z_s;
370 N = -1 * (Dn) * dx * z_n + vn * (1 - rn) * dx * z_n;
389 C = (Dw + vw * rw) * dy * z_w +
390 (De + ve * re) * dy * z_e +
391 (Ds + vs * rs) * dx * z_s +
392 (Dn + vn * rn) * dx * z_n + Az * z * R / data->
dt - q / nf;
395 V = (cs + cg_start * Az * z * R / data->
dt + q / nf * cin);
410 G_debug(6,
"N_callback_solute_transport_2d: called [%i][%i]", row, col);
626 "N_calc_solute_transport_transmission_2d: calculating transmission boundary");
628 for (j = 0; j < rows; j++) {
629 for (i = 0; i <
cols; i++) {
663 c = c / (double)count;
665 if (c > 0 || c == 0 || c < 0)
696 double disp_xx, disp_yy, disp_xy;
704 "N_calc_solute_transport_disptensor_2d: calculating the dispersivity tensor");
706 for (j = 0; j < rows; j++) {
707 for (i = 0; i <
cols; i++) {
715 vx = (grad.
WC + grad.
EC) / 2;
716 vy = (grad.
NC + grad.
SC) / 2;
717 vv = sqrt(vx * vx + vy * vy);
720 disp_xx = data->
al * vx * vx / vv + data->
at * vy * vy / vv;
721 disp_yy = data->
at * vx * vx / vv + data->
al * vy * vy / vv;
722 disp_xy = (data->
al - data->
at) * vx * vy / vv;
726 "N_calc_solute_transport_disptensor_2d: [%i][%i] disp_xx %g disp_yy %g disp_xy %g",
727 i, j, disp_xx, disp_yy, disp_xy);
755 int cols, rows, depths;
757 double vx, vy, vz, vv;
759 double disp_xx, disp_yy, disp_zz, disp_xy, disp_xz, disp_yz;
768 "N_calc_solute_transport_disptensor_3d: calculating the dispersivity tensor");
770 for (k = 0; k < depths; k++) {
771 for (j = 0; j < rows; j++) {
772 for (i = 0; i <
cols; i++) {
782 vx = (grad.
WC + grad.
EC) / 2;
783 vy = (grad.
NC + grad.
SC) / 2;
784 vz = (grad.
BC + grad.
TC) / 2;
785 vv = sqrt(vx * vx + vy * vy + vz * vz);
789 data->
al * vx * vx / vv + data->
at * vy * vy / vv +
790 data->
at * vz * vz / vv;
792 data->
at * vx * vx / vv + data->
al * vy * vy / vv +
793 data->
at * vz * vz / vv;
795 data->
at * vx * vx / vv + data->
at * vy * vy / vv +
796 data->
al * vz * vz / vv;
797 disp_xy = (data->
al - data->
at) * vx * vy / vv;
798 disp_xz = (data->
al - data->
at) * vx * vz / vv;
799 disp_yz = (data->
al - data->
at) * vy * vz / vv;
803 "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 ",
804 i, j, k, disp_xx, disp_yy, disp_zz, disp_xy, disp_xz,
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
N_solute_transport_data3d * N_alloc_solute_transport_data3d(int cols, int rows, int depths)
Alllocate memory for the solute transport data structure in three dimensions.
void G_free(void *buf)
Free allocated memory.
Matrix entries for a mass balance 5/7/9 star system.
double N_calc_geom_mean(double a, double b)
Calculate the geometrical mean of values a and b.
#define N_CELL_TRANSMISSION
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Gradient between the cells in X and Y direction.
void N_calc_solute_transport_transmission_2d(N_solute_transport_data2d *data)
Compute the transmission boundary condition in 2d.
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.
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
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.
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...
N_gradient_field_3d * grad
Geometric information about the structured grid.
N_gradient_field_2d * N_alloc_gradient_field_2d(int cols, int rows)
Allocate a N_gradient_field_2d.
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
void N_calc_solute_transport_disptensor_3d(N_solute_transport_data3d *data)
Compute the dispersivity tensor based on the solute transport data in 3d.
double N_exp_upwinding(double sprod, double distance, double D)
exponential upwinding stabilization algorithm
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
void N_calc_solute_transport_disptensor_2d(N_solute_transport_data2d *data)
Compute the dispersivity tensor based on the solute transport data in 2d.
void N_free_solute_transport_data2d(N_solute_transport_data2d *data)
Release the memory of the solute transport data structure in two dimensions.
Gradient between the cells in X, Y and Z direction.
void N_free_solute_transport_data3d(N_solute_transport_data3d *data)
Release the memory of the solute transport data structure in three dimensions.
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.
N_gradient_field_2d * grad
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 postion col, row is of type null, otherwise 0...
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.
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.
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]...
int G_debug(int level, const char *msg,...)
Print debugging message.
N_gradient_field_3d * N_alloc_gradient_field_3d(int cols, int rows, int depths)
Allocate a N_gradient_field_3d.
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.
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.
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
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.
N_solute_transport_data2d * N_alloc_solute_transport_data2d(int cols, int rows)
Alllocate memory for the solute transport data structure in two dimensions.
double N_full_upwinding(double sprod, double distance, double D)
full upwinding stabilization algorithm