2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../classes/Node.h"
5 #include "../shared/shared.h"
6 #include "../modules/modules.h"
7 #include "../classes/Inputs2/TransientInput2.h"
16 int sedimentlimit_flag;
24 bool isefficientlayer;
29 char** requestedoutputs = NULL;
34 iomodel->
FindConstant(&hydrology_model,
"md.hydrology.model");
39 iomodel->
FetchData(&sedimentlimit_flag,
"md.hydrology.sedimentlimit_flag" );
40 iomodel->
FetchData(&transfer_flag,
"md.hydrology.transfer_flag" );
41 iomodel->
FetchData(&unconfined_flag,
"md.hydrology.unconfined_flag" );
42 iomodel->
FetchData(&penalty_lock,
"md.hydrology.penalty_lock" );
43 iomodel->
FetchData(&hydro_maxiter,
"md.hydrology.max_iter" );
44 iomodel->
FetchData(&hydroslices,
"md.hydrology.steps_per_step");
45 iomodel->
FetchData(&averaging_method,
"md.hydrology.averaging");
46 iomodel->
FetchData(&isefficientlayer,
"md.hydrology.isefficientlayer");
47 iomodel->
FetchData(&penalty_factor,
"md.hydrology.penalty_factor" );
48 iomodel->
FetchData(&rel_tol,
"md.hydrology.rel_tol" );
63 iomodel->
FetchData(&leakagefactor,
"md.hydrology.leakage_factor");
66 if(sedimentlimit_flag==1){
67 iomodel->
FetchData(&sedimentlimit,
"md.hydrology.sedimentlimit");
76 iomodel->
FindConstant(&requestedoutputs,&numoutputs,
"md.hydrology.requested_outputs");
79 iomodel->
DeleteData(&requestedoutputs,numoutputs,
"md.hydrology.requested_outputs");
83 bool isefficientlayer;
87 iomodel->
FindConstant(&hydrology_model,
"md.hydrology.model");
93 iomodel->
FindConstant(&isefficientlayer,
"md.hydrology.isefficientlayer");
100 element->
Update(inputs2,i,iomodel,analysis_counter,analysis_type,
P1Enum);
116 if(isefficientlayer){
126 iomodel->
FindConstant(&hydrology_model,
"md.hydrology.model");
132 iomodel->
FetchData(2,
"md.mesh.vertexonbase",
"md.mesh.vertexonsurface");
135 iomodel->
DeleteData(2,
"md.mesh.vertexonbase",
"md.mesh.vertexonsurface");
141 iomodel->
FindConstant(&hydrology_model,
"md.hydrology.model");
150 iomodel->
FindConstant(&hydrology_model,
"md.hydrology.model");
154 iomodel->
FetchData(1,
"md.mesh.vertexonbase");
166 else if(reCast<int>(iomodel->
Data(
"md.mesh.vertexonbase")[i])){
173 iomodel->
DeleteData(1,
"md.mesh.vertexonbase");
198 basalelement = element;
201 if(!element->
IsOnBase())
return NULL;
210 if(!thawed_element) {
219 bool active_element,isefficientlayer;
231 IssmDouble* basis = xNew<IssmDouble>(numnodes);
244 if(isefficientlayer){
251 for(
int ig=gauss -> begin();ig<gauss->
end();ig++){
252 gauss -> GaussPoint(ig);
253 basalelement -> JacobianDeterminant(&Jdet,xyz_list,gauss);
255 sediment_transmitivity =
SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
256 sediment_storing =
SedimentStoring(basalelement,gauss,sed_head_input,base_input);
259 D_scalar=sediment_transmitivity*gauss->
weight*Jdet;
261 if(dt!=0.) D_scalar=D_scalar*dt;
264 GetB(B,basalelement,xyz_list,gauss);
273 D_scalar=sediment_storing*gauss->
weight*Jdet;
281 if(isefficientlayer){
285 D_scalar=dt*transfer*gauss->
weight*Jdet;
296 xDelete<IssmDouble>(xyz_list);
297 xDelete<IssmDouble>(B);
298 xDelete<IssmDouble>(basis);
317 basalelement = element;
320 if(!element->
IsOnBase())
return NULL;
329 if(!thawed_element) {
338 bool active_element,isefficientlayer;
341 int hydrologysubstepping;
349 Input2 *active_element_input = NULL;
350 Input2 *old_wh_input = NULL;
351 Input2 *dummy_input = NULL;
352 Input2 *surface_runoff_input = NULL;
359 IssmDouble* basis = xNew<IssmDouble>(numnodes);
384 if(smbsubstepping==1){
388 else if(smbsubstepping>1 && smbsubstepping<=hydrologysubstepping){
396 surface_runoff_input=xDynamicCast<Input2*>(dummy_input);
_assert_(surface_runoff_input);
400 if(isefficientlayer){
410 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
414 sediment_transmitivity =
SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
417 if(!isefficientlayer){
419 if(surface_runoff_input) surface_runoff_input->
GetInputValue(&runoff_value,gauss);
420 else runoff_value = 0.;
421 scalar = Jdet*gauss->
weight*(water_load+runoff_value);
423 if(dt!=0.) scalar = scalar*dt;
424 for(
int i=0;i<numnodes;i++){
425 pe->
values[i]+=scalar*basis[i];
432 if(surface_runoff_input)surface_runoff_input->
GetInputValue(&runoff_value,gauss);
433 else runoff_value = 0.;
434 scalar = Jdet*gauss->
weight*(water_load+runoff_value);
436 if(dt!=0.) scalar = scalar*dt;
437 for(
int i=0;i<numnodes;i++){
438 pe->
values[i]+=scalar*basis[i];
446 sediment_storing =
SedimentStoring(basalelement,gauss,sed_head_input,base_input);
447 if(isefficientlayer){
455 scalar = Jdet*gauss->
weight*((water_head*sediment_storing)+(dt*transfer));
457 for(
int i=0;i<numnodes;i++)pe->
values[i]+=scalar*basis[i];
460 scalar = Jdet*gauss->
weight*(water_head*sediment_storing);
462 for(
int i=0;i<numnodes;i++)pe->
values[i]+=scalar*basis[i];
467 xDelete<IssmDouble>(xyz_list);
468 xDelete<IssmDouble>(basis);
491 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
495 for(
int i=0;i<numnodes;i++){
496 B[numnodes*0+i] = dbasis[0*numnodes+i];
497 B[numnodes*1+i] = dbasis[1*numnodes+i];
501 xDelete<IssmDouble>(dbasis);
507 _error_(
"Not implemented yet");
521 basalelement = element;
534 IssmDouble* values = xNew<IssmDouble>(numnodes);
535 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
536 IssmDouble* residual = xNew<IssmDouble>(numnodes);
538 for(
int i=0;i<numnodes;i++){
539 values[i] =solution[doflist[i]];
540 if(xIsNan<IssmDouble>(values[i]))
_error_(
"NaN found in solution vector");
541 if(xIsInf<IssmDouble>(values[i]))
_error_(
"Inf found in solution vector");
551 IssmDouble* thickness = xNew<IssmDouble>(numnodes);
552 IssmDouble* base = xNew<IssmDouble>(numnodes);
553 IssmDouble* transmitivity = xNew<IssmDouble>(numnodes);
565 kappa=kmax*pow(10.,penalty_factor);
567 for(
int i=0;i<numnodes;i++){
569 if(values[i]>h_max) {
570 residual[i] = kappa*(values[i]-h_max);
576 pressure[i]=(rho_ice*g*thickness[i])-(rho_freshwater*g*(values[i]-base[i]));
578 xDelete<IssmDouble>(thickness);
579 xDelete<IssmDouble>(base);
580 xDelete<IssmDouble>(transmitivity);
589 xDelete<IssmDouble>(values);
590 xDelete<IssmDouble>(residual);
591 xDelete<IssmDouble>(pressure);
592 xDelete<int>(doflist);
608 IssmDouble base_elev,prestep_head,water_sheet;
616 switch(unconf_scheme){
618 sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
623 water_sheet=
max(0.0,(prestep_head-(base_elev-sediment_thickness)));
631 storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
634 yield=sediment_porosity;
635 sediment_storing=yield+(storing-yield)/(1+exp(-2*expfac*(water_sheet-0.99*sediment_thickness)));
638 _error_(
"UnconfinedFlag is 0 or 1");
640 return sediment_storing;
649 IssmDouble base_elev,prestep_head,water_sheet;
653 SedTrans_input->
GetInputValue(&FullLayer_transmitivity,gauss);
655 switch(unconf_scheme){
657 sediment_transmitivity=FullLayer_transmitivity;
662 water_sheet=
max(0.0,(prestep_head-(base_elev-sediment_thickness)));
665 sediment_transmitivity=FullLayer_transmitivity*
min(water_sheet/sediment_thickness,1.);
666 if (sediment_transmitivity==0){
667 sediment_transmitivity=1.0e-20;
672 _error_(
"UnconfinedFlag is 0 or 1");
674 return sediment_transmitivity;
699 h_max=((rho_ice*thickness)/rho_water)+bed;
702 _error_(
"Using normal stress not supported yet");
705 _error_(
"no case higher than 3 for SedimentlimitFlag");
718 switch(transfermethod){
728 _error_(
"no case higher than 1 for the Transfer method");
740 switch(transfermethod){
749 transfer=+epl_head*leakage;
752 _error_(
"no case higher than 1 for the Transfer method");
766 element_active =
true;
769 element_active =
false;
783 basalelement = element;
793 IssmDouble* meltingrate = xNew<IssmDouble>(numnodes);
794 IssmDouble* groundedice = xNew<IssmDouble>(numnodes);
800 for(
int i=0;i<numnodes;i++){
801 if ((meltingrate[i]<=0.0) && (groundedice[i]>0)){
813 xDelete<IssmDouble>(meltingrate);
814 xDelete<IssmDouble>(groundedice);
826 element_active =
true;
829 element_active =
false;
843 basalelement = element;
854 IssmDouble* active = xNew<IssmDouble>(numnodes);
861 for(
int i=0;i<numnodes;i++) flag+=active[i];
865 for(
int i=0;i<numnodes;i++){
870 else if(active_element){
871 for(
int i=0;i<numnodes;i++){
879 xDelete<IssmDouble>(active);