Changeset 18745

Timestamp:

11/06/14 10:02:00 (10 years ago)

Author:

bdef

Message:

CHG: Now applying exponential heavyside approximation to the two forks of transfer term computation

Location:

issm/trunk-jpl/src/c/analyses

Files:

• HydrologyDCEfficientAnalysis.cpp (modified) (8 diffs)
• HydrologyDCInefficientAnalysis.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

@@ -426 +426 @@
         IssmDouble epl_head,sediment_head;
         IssmDouble leakage,transfer;
+        IssmDouble continuum, factor;
         HydrologyDCInefficientAnalysis* inefanalysis = NULL;
 
@@ -448 +449 @@
                 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
                 
-                if(epl_head>sediment_head){  
-                        /* if(sediment_head>=hmax){ */
-                        /*      transfer=0.0; */
-                        /* } */
-                        /* else{ */
-                        /*      transfer=(leakage); */
-                        /* } */
-                        transfer=leakage*(1/(1+exp(-20.0*(hmax-sediment_head))));
-                }
-                else{
-                        transfer=(leakage);
-                        
-                }
+                //Computing continuum function to apply to transfer term, transfer is null only if
+                // epl_head>sediment_head AND sediment_head>h_max
+                continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head))));
+                factor=max(continuum,1.0);
+                transfer=leakage*factor;
                 break;
         default:
@@ -474 +467 @@
         IssmDouble epl_head,sediment_head;
         IssmDouble leakage,transfer;
+        IssmDouble continuum, factor;
 
         HydrologyDCInefficientAnalysis* inefanalysis = NULL;
@@ -497 +491 @@
                 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
 
-                if(epl_head>sediment_head){  
-                        /* if(sediment_head>=hmax){ */
-                        /*      transfer=0.0; */
-                        /* } */
-                        /* else{ */
-                        /*      transfer=(sediment_head*leakage); */
-                        /* } */
-                        transfer=sediment_head*leakage*(1/(1+exp(-20.0*(hmax-sediment_head))));
-                }
-                else{
-                        transfer=(sediment_head*leakage);
-                }
+                //Computing continuum function to apply to transfer term, transfer is null only if
+                // epl_head>sediment_head AND sediment_head>h_max
+                continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head))));
+                factor=max(continuum,1.0);
+                transfer=sediment_head*leakage*factor;
+
                 break;
         default:
@@ -738 +726 @@
 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
         /*Constants*/
-
         int      domaintype;
         Element*   basalelement=NULL;
@@ -760 +747 @@
                 
         basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
-
         bool active_element;
-        
         Input*  active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);           
         active_element_input->GetInputValue(&active_element);
-
         
         for(int i=0;i<numnodes;i++) flag+=active[i];
@@ -775 +759 @@
         }
         else if(active_element){
-                /*Checking Stuff*/
                 for(int i=0;i<numnodes;i++){
                         active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
@@ -781 +764 @@
         }
         else{
-                /*Do not do anything: at least one node was active for this element but this element is not solved for*/
+                /*Do not do anything: at least one node is active for this element but this element is not solved for*/
         }
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

@@ -593 +593 @@
         IssmDouble epl_head,sediment_head;
         IssmDouble leakage,transfer;
-
+        IssmDouble continuum, factor;
+        
         element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
 
@@ -611 +612 @@
 
                 hmax=GetHydrologyDCInefficientHmax(element, gauss, thick_input, base_input);
-                
-                if(epl_head>sediment_head){  
-                        /* if(sediment_head>=hmax){ */
-                        /*      transfer=0.0; */
-                        /* } */
-                        /* else{ */
-                        /*      transfer=(leakage); */
-                        /* } */
-                        transfer=leakage*(1/(1+exp(-20.0*(hmax-sediment_head))));
-                }
-                else{
-                        transfer=(leakage);
-                }               
+        
+                //Computing continuum function to apply to transfer term, transfer is null only if
+                //epl_head>sediment_head AND sediment_head>h_max
+                continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head))));
+                factor=max(continuum,1.0);
+                transfer=leakage*factor;
+
                 break;
         default:
@@ -637 +632 @@
         IssmDouble epl_head,sediment_head;
         IssmDouble leakage,transfer;
-
+        IssmDouble continuum, factor;
+        
         element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
 
@@ -656 +652 @@
                 hmax=GetHydrologyDCInefficientHmax(element, gauss, thick_input, base_input);
                 
-                if(epl_head>sediment_head){  
-                        /* if(sediment_head>=hmax){ */
-                        /*      transfer=0.0; */
-                        /* } */
-                        /* else{ */
-                        /*      transfer=(epl_head*leakage); */
-                        /* } */
-                        transfer=epl_head*leakage*(1/(1+exp(-20.0*(hmax-sediment_head))));
-                }
-                else{
-                        transfer=(epl_head*leakage);
-                }
+                //Computing continuum function to apply to transfer term, transfer is null only if
+                //epl_head>sediment_head AND sediment_head>h_max
+                continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head))));
+                factor=max(continuum,1.0);
+                transfer=epl_head*leakage*factor;
                 break;
         default: