00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "richards.h"
00013 #include "globalt.h"
00014
00015
00016 richards::richards()
00017 {
00018
00019 alpha=0.0;
00020 n=0.0;
00021 m=0.0;
00022
00023 dim=0;
00024
00025
00026 kksxx=0.0;
00027 kksxy=0.0;
00028 kksxz=0.0;
00029 kksyy=0.0;
00030 kksyz=0.0;
00031 kkszz=0.0;
00032
00033
00034 thetas=0.0;
00035
00036 thetar=0.0;
00037
00038 }
00039
00040 richards::~richards()
00041 {
00042
00043 }
00044
00045
00046
00047
00048
00049
00050
00051
00052 void richards::read(XFILE *in)
00053 {
00054
00055
00056 dim = Tt->give_dimension (0);
00057 if (dim<1 || dim>3){
00058 print_err("wrong dimension of problem solved (%ld)",__FILE__,__LINE__,__func__,dim);
00059 }
00060
00061 switch (dim) {
00062 case 1:
00063 xfscanf (in,"%lf",&kksxx);
00064 break;
00065 case 2:
00066 xfscanf (in,"%lf %lf %lf",&kksxx,&kksxy,&kksyy);
00067 break;
00068 case 3:
00069 xfscanf (in,"%lf %lf %lf %lf %lf %lf",&kksxx,&kksxy,&kksxz,&kksyy,&kksyz,&kkszz);
00070 break;
00071 }
00072
00073
00074
00075
00076
00077 xfscanf (in,"%lf %lf %lf %lf %lf %lf",&alpha,&n,&m,&thetas,&thetar,&storage);
00078 }
00079
00080
00081
00082
00083
00084
00085
00086
00087 void richards::print(FILE *out)
00088 {
00089 switch (dim) {
00090 case 1:
00091 fprintf (out,"%le",kksxx);
00092 break;
00093 case 2:
00094 fprintf (out,"%le %le %le",kksxx,kksxy,kksyy);
00095 break;
00096 case 3:
00097 fprintf (out,"%le %le %le %le %le %le",kksxx,kksxy,kksxz,kksyy,kksyz,kkszz);
00098 break;
00099 }
00100
00101
00102
00103
00104
00105 fprintf (out," %le %le %le %le %le\n",alpha,n,m,thetas,thetar);
00106 }
00107
00108
00109
00110
00111
00112
00113
00114
00115 double richards::kk_value (double h)
00116 {
00117 if (h < 0){
00118 return (pow(1 - pow(-(alpha*h),m*n)/pow(1 + pow(-(alpha*h),n),m),2)/pow(1 + pow(-(alpha*h),n),m/2.0));
00119 }
00120 else{
00121 return 1.0;
00122 }
00123
00124 }
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134 double richards::dkkdh_value (double h)
00135 {
00136
00137 if (h < 0){
00138 return ((alpha*pow(-(alpha*h),-1 + n)*pow(1 + pow(-(alpha*h),n),-1 - m/2.)*pow(1 - pow(-(alpha*h),m*n)/pow(1 + pow(-(alpha*h),n),m),2)*m*n)/2. + (2*(1 - pow(-(alpha*h),m*n)/pow(1 + pow(-(alpha*h),n),m))*(-(alpha*pow(-(alpha*h),-1 + n + m*n)*pow(1 + pow(-(alpha*h),n),-1 - m)*m*n) + (alpha*pow(-(alpha*h),-1 + m*n)*m*n)/pow(1 + pow(-(alpha*h),n),m)))/pow(1 + pow(-(alpha*h),n),m/2.));
00139 }
00140 else{
00141 return 0.0;
00142 }
00143
00144 }
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206 double richards::theta_val (double h)
00207 {
00208 if (h < 0){
00209 return (thetas - thetar)/pow(1 + pow(-(alpha*h),n),m) + thetar ;
00210 }
00211 else{
00212 return thetas;
00213 }
00214 }
00215
00216
00217
00218
00219
00220
00221
00222
00223 vector richards::darcian_flux (long ipp)
00224 {
00225 double h;
00226 vector gradH(dim);
00227 vector flux(dim);
00228 matrix K(dim,dim);
00229 int i, j;
00230 double Kr;
00231
00232 h = Tm->ip[ipp].av[0] ;
00233
00234 Kr = kk_value(h) ;
00235
00236 switch (dim) {
00237 case 1:
00238 K(0,0) = kksxx;
00239 break ;
00240 case 2:
00241 K(0,0) = kksxx;
00242 K(0,1) = kksxy;
00243 K(1,0) = kksxy;
00244 K(1,1) = kksyy;
00245 break;
00246 case 3:
00247 K(0,0) = kksxx;
00248 K(0,1) = kksxy;
00249 K(0,2) = kksxz;
00250 K(1,0) = kksxy;
00251 K(1,1) = kksyy;
00252 K(1,2) = kksyz;
00253 K(2,0) = kksxz;
00254 K(2,1) = kksyz;
00255 K(2,2) = kkszz;
00256 break;
00257 }
00258
00259
00260 for (i = 1; i <= dim; i++) {
00261 for (j = 1; j <= dim; j++) {
00262 K(i-1,j-1) = K(i-1,j-1)*Kr*(-1.0) ;
00263 }
00264 }
00265
00266 for (i = 1; i <= dim; i++) {
00267 gradH(i-1) = Tm->ip[ipp].grad[0][i-1];
00268 }
00269
00270 mxv(K, gradH, flux) ;
00271
00272 return flux ;
00273
00274
00275 }
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 double richards::c_value (double h)
00288 {
00289 if (h < 0){
00290 return alpha*pow(-(alpha*h),-1 + n)*pow(1 + pow(-(alpha*h),n),-1 - m)*m*n*(thetas - thetar) +
00291 storage*theta_val(h)/thetas;
00292 }
00293 else{
00294 return storage;
00295 }
00296 }
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 double richards::dcdh_value (double h)
00308 {
00309 if (h < 0){
00310 return -(pow(alpha,2)*pow(-(alpha*h),-2 + n)*
00311 pow(1 + pow(-(alpha*h),n),-1 - m)*m*(-1 + n)*n*
00312 (thetas - thetar)) - pow(alpha,2)*pow(-(alpha*h),-2 + 2*n)*
00313 pow(1 + pow(-(alpha*h),n),-2 - m)*(-1 - m)*m*pow(n,2)*
00314 (thetas - thetar) + (alpha*pow(-(alpha*h),-1 + n)*
00315 pow(1 + pow(-(alpha*h),n),-1 - m)*m*n*storage*
00316 (thetas - thetar))/thetas;
00317 }
00318 else{
00319 return 0.0;
00320 }
00321 }
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334 void richards::matcond (matrix &d,long ri,long ci,long ipp)
00335 {
00336 long n;
00337 n = d.n;
00338
00339 switch (n){
00340 case 1:{
00341 matcond1d (d,ri,ci,ipp);
00342 break;
00343 }
00344 case 2:{
00345 matcond2d (d,ri,ci,ipp);
00346 break;
00347 }
00348 case 3:{
00349 matcond3d (d,ri,ci,ipp);
00350 break;
00351 }
00352 default:{
00353 print_err("unknown number of flux components is required",__FILE__,__LINE__,__func__);
00354 }
00355 }
00356 }
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 void richards::matcond1d (matrix &d,long ri,long ci,long ipp)
00368 {
00369 double h,k;
00370
00371
00372 h = Tm->ip[ipp].av[0];
00373
00374
00375 k = kk_value (h);
00376
00377 d[0][0] = k*kksxx;
00378 }
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 void richards::matcond2d (matrix &d,long ri,long ci,long ipp)
00389 {
00390 double h,k;
00391
00392
00393 h = Tm->ip[ipp].av[0];
00394
00395
00396 k = kk_value (h);
00397
00398 fillm(0.0,d);
00399
00400 d[0][0] = k*kksxx; d[0][1] = k*kksxy;
00401 d[1][0] = k*kksxy; d[1][1] = k*kksyy;
00402
00403
00404 }
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415 void richards::matcond3d (matrix &d,long ri,long ci,long ipp)
00416 {
00417 double h,k;
00418
00419
00420 h = Tm->ip[ipp].av[0];
00421
00422
00423 k = kk_value (h);
00424
00425 fillm(0.0,d);
00426
00427 d[0][0] = k*kksxx; d[0][1] = k*kksxy; d[0][2]=k*kksxz;
00428 d[1][0] = k*kksxy; d[1][1] = k*kksyy; d[1][2]=k*kksyz;
00429 d[2][0] = k*kksxz; d[2][1] = k*kksyz; d[2][2]=k*kkszz;
00430 }
00431
00432 void richards::matcond2 (matrix &d,long ri,long ci,long ipp)
00433 {
00434 long n;
00435 n = d.n;
00436
00437 switch (n){
00438 case 1:{
00439 matcond1d_2 (d,ri,ci,ipp);
00440 break;
00441 }
00442 case 2:{
00443 matcond2d_2 (d,ri,ci,ipp);
00444 break;
00445 }
00446 case 3:{
00447 matcond3d_2 (d,ri,ci,ipp);
00448 break;
00449 }
00450 default:{
00451 print_err("unknown number of components is required",__FILE__,__LINE__,__func__);
00452 }
00453 }
00454 }
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466 void richards::matcond1d_2 (matrix &d,long ri,long ci,long ipp)
00467 {
00468 double h;
00469
00470
00471 h = Tm->ip[ipp].av[0];
00472
00473 fillm(0.0,d);
00474 d[0][0]=0.0-kksxx*dkkdh_value (h);
00475 }
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489 void richards::matcond2d_2 (matrix &d,long ri,long ci,long ipp)
00490 {
00491 double h;
00492
00493
00494 h = Tm->ip[ipp].av[0];
00495
00496 fillm(0.0,d);
00497 d[0][0]=0.0-kksxy*dkkdh_value (h);
00498 d[0][1]=0.0-kksyy*dkkdh_value (h);
00499
00500
00501 }
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513 void richards::matcond3d_2 (matrix &d,long ri,long ci,long ipp)
00514 {
00515 double h;
00516
00517
00518 h = Tm->ip[ipp].av[0];
00519
00520 fillm(0.0,d);
00521 d[0][0]=0.0 - kksxz*dkkdh_value (h);
00522 d[0][1]=0.0 - kksyz*dkkdh_value (h);
00523 d[0][2]=0.0 - kkszz*dkkdh_value (h);
00524 }
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535 void richards::matcap (double &c,long ri,long ci,long ipp)
00536 {
00537 double h;
00538
00539
00540 h = Tm->ip[ipp].av[0];
00541
00542
00543 c = c_value (h);
00544
00545 }
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559 void richards::give_dof_names(namevart *dofname, long ntm)
00560 {
00561 if (ntm < 1)
00562 {
00563 print_err("the model defines %ld unknowns but number of transported media is %ld",
00564 __FILE__, __LINE__, __func__, 1, ntm);
00565 abort();
00566 }
00567 dofname[0] = trf_hydraulic_head;
00568 }
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579 double richards::taylor_error(long ipp)
00580 {
00581 double hn;
00582 double hp;
00583
00584
00585 hn = Tm->ip[ipp].av[0];
00586
00587
00588 hp = Tm->ip[ipp].pv[0];
00589
00590 return theta_val(hp) + c_value(hp)*(hp-hn) - theta_val(hn) ;
00591
00592 }
00593
00594
00595
00596
00597