Changeset 27518


Ignore:
Timestamp:
01/19/23 03:50:29 (2 years ago)
Author:
bdef
Message:

BUG: fix to the split between layers of moulin input for DoCo

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r27177 r27518  
    297297        IssmDouble dt,scalar,water_head;
    298298        IssmDouble water_load,transfer,runoff_value;
    299         IssmDouble epl_storing,epl_transmitivity;
     299        IssmDouble epl_storing;  //,epl_transmitivity;
    300300        IssmDouble Jdet,time;
    301301        IssmDouble residual,connectivity;
     302        IssmDouble active_node;
    302303
    303304        IssmDouble *xyz_list             = NULL;
    … …  
    313314        ElementVector* pe    = basalelement->NewElementVector();
    314315        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     316
    315317
    316318        /*Retrieve all inputs and parameters*/
    … …  
    353355                basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss);
    354356                basalelement ->NodalFunctions(basis,gauss);
     357                //epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input);
    355358                epl_storing     = EplStoring(basalelement,gauss,epl_thick_input);
    356                 epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input);
    357359
    358360                /*Loading term*/
    … …  
    362364                scalar = Jdet*gauss->weight*(water_load+runoff_value);
    363365                if(dt!=0.) scalar = scalar*dt;
    364                 for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    365 
     366                for(int i=0;i<numnodes;i++){
     367                        //This is the original
     368                        pe->values[i]+=scalar*basis[i];
     369                        //This is the noded version
     370                        /* basalelement->GetInputValue(&active_node,basalelement->nodes[i],HydrologydcMaskEplactiveNodeEnum); */
     371                        /* if(!reCast<bool>(active_node)){ */
     372                        /*      pe->values[i]+=scalar*basis[i]; */
     373                        //}
     374                        //if(basalelement->nodes[i]->Sid()==42)_printf_("EPL uni Input "<<scalar*basis[i]<<"\n");
     375                }
    366376                /*Transient and transfer terms*/
    367377                if(dt!=0.){
    … …  
    379389        for(int iv=0;iv<numvertices;iv++){
    380390                gauss->GaussVertex(iv);
    381                 epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input);
     391                //epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input);
    382392                connectivity = IssmDouble(basalelement->VertexConnectivity(iv));
    383393                residual_input->GetInputValue(&residual,gauss);

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

    r27177 r27518  
    326326        /*Intermediaries */
    327327        bool       active_element,isefficientlayer;
    328         int        smb_model,smbsubstepping;
    329         int        hydrologysubstepping,smb_averaging;
    330         IssmDouble dt,scalar,sediment_storing;
    331         IssmDouble water_head,sediment_transmitivity;
     328        int        smb_model,smb_averaging;
     329        int        smbsubstepping, hydrologysubstepping;
     330        IssmDouble dt,scalar,water_head;
     331        IssmDouble sediment_storing,sediment_transmitivity;
    332332        IssmDouble water_load,runoff_value,transfer;
    333333        IssmDouble Jdet,time;
     334        IssmDouble active_node;
    334335
    335336        IssmDouble *xyz_list             = NULL;
    … …  
    401402                        scalar = Jdet*gauss->weight*(water_load+runoff_value);
    402403                        if(dt!=0.) scalar = scalar*dt;
    403                         for(int i=0;i<numnodes;i++){
    404                                 pe->values[i]+=scalar*basis[i];
    405                         }
     404                        for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    406405                }
    407406                else{
    … …  
    414413                                if(dt!=0.) scalar = scalar*dt;
    415414                                for(int i=0;i<numnodes;i++){
     415                                        //This is the original
    416416                                        pe->values[i]+=scalar*basis[i];
     417                                        //This is the noded version
     418                                        /* basalelement->GetInputValue(&active_node,basalelement->nodes[i],HydrologydcMaskEplactiveNodeEnum); */
     419                                        /* if(!reCast<bool>(active_node)){ */
     420                                        /*      pe->values[i]+=scalar*basis[i]; */
     421                                        /*      //if(basalelement->nodes[i]->Sid()==42)_printf_("IDS uni Input "<<scalar*basis[i]<<"\n"); */
     422                                        //}
    417423                                }
    418424                        }

  • issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp

    r26144 r27518  
    104104/*}}}*/
    105105void    Moulin::Echo(void){/*{{{*/
    106         this->DeepEcho();
     106
     107        _printf_("Moulin:\n");
     108        _printf_("   id: " << id << "\n");
     109        hnode->Echo();
     110        hvertex->Echo();
     111        helement->Echo();
     112        _printf_("   parameters\n");
     113        parameters->Echo();
     114        //this->DeepEcho();
    107115}
    108116/*}}}*/
    … …  
    143151void  Moulin::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
    144152
    145         /*Take care of hooking up all objects for this load, ie links the objects in the hooks to their respective 
     153        /*Take care of hooking up all objects for this load, ie links the objects in the hooks to their respective
    146154         * datasets, using internal ids and offsets hidden in hooks: */
    147155        hnode->configure(nodesin);
    … …  
    204212                        break;
    205213                case HydrologyDCInefficientAnalysisEnum:
    206                         pe = CreatePVectorHydrologyDCInefficient();
     214                        pe = this->CreatePVectorHydrologyDCInefficient();
    207215                        break;
    208216                case HydrologyDCEfficientAnalysisEnum:
    209                         pe = CreatePVectorHydrologyDCEfficient();
     217                        pe = this->CreatePVectorHydrologyDCEfficient();
    210218                        break;
    211219                default:
    … …  
    427435         * mesh), don't add the moulin input a second time*/
    428436        if(node->IsClone()) return NULL;
    429         bool isefficientlayer;
     437        bool isefficientlayer, active_element;
    430438        IssmDouble moulin_load,dt;
    431439        IssmDouble epl_active;
    432 
    433440        /*Initialize Element matrix*/
    434441        ElementVector* pe=new ElementVector(&node,1,this->parameters);
    … …  
    437444        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    438445        parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    439         // Test version input in EPL when active
     446
     447        //Test version input in EPL when active
    440448        if(isefficientlayer){
    441                 this->element->GetInputValue(&epl_active,node,HydrologydcMaskEplactiveNodeEnum);
    442                 if(reCast<bool>(epl_active)){
    443                         pe->values[0]=moulin_load*0.0;
     449                this->element->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum);
     450                if(!active_element){
     451                        /* this->element->GetInputValue(&epl_active,node,HydrologydcMaskEplactiveNodeEnum); */
     452                        /* if(reCast<bool>(epl_active))pe->values[0]=0.0; */
     453                        /* else { */
     454                        pe->values[0]=moulin_load*dt;
     455                        /*      if (moulin_load>0)_printf_("MoulinInput in Sed is "<<pe->values[0]<<"\n"); */
     456                        /*      if (moulin_load>0)pe->Echo(); */
     457                        /* } */
     458                        //if (node->Sid()==4)_printf_("MoulinInput in Sed is "<<moulin_load*dt<<"\n");
    444459                }
    445                 else{
    446                         pe->values[0]=moulin_load*dt;
    447                 }
    448         }
    449         else{
    450                 pe->values[0]=moulin_load*dt;
    451         }
     460                else pe->values[0]=0.0;
     461        }
     462        else pe->values[0]=moulin_load*dt;
     463
     464        //Test only input in sed
     465        /* pe->values[0]=moulin_load*dt; */
     466
    452467        /*Clean up and return*/
    453468        return pe;
    … …  
    458473        /*If this node is not the master node (belongs to another partition of the
    459474         * mesh), don't add the moulin input a second time*/
     475
    460476        if(node->IsClone()) return NULL;
    461         if(!this->node->IsActive()) return NULL;
    462         IssmDouble moulin_load,dt;
    463477        ElementVector* pe=new ElementVector(&node,1,this->parameters);
    464478
    465         this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
    466         parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    467 
    468         pe->values[0]=moulin_load*dt;
    469         /*Clean up and return*/
     479        //Test Input in epl if active
     480        /* IssmDouble epl_active; */
     481        /* this->element->GetInputValue(&epl_active,node,HydrologydcMaskEplactiveNodeEnum); */
     482        /* //if(node->Sid()==4)_printf_("Activity is "<<epl_active<<" \n"); */
     483        /* if(reCast<bool>(epl_active)){ */
     484        /*      IssmDouble moulin_load,dt; */
     485        /*      this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum); */
     486        /*      parameters->FindParam(&dt,TimesteppingTimeStepEnum); */
     487        /*      pe->values[0]=moulin_load*dt; */
     488        /*      if (moulin_load>0)_printf_("MoulinInput in Epl is "<<pe->values[1]<<"\n"); */
     489
     490        /* } */
     491        /*      else{ */
     492        /*              pe->values[0]=0.0; */
     493        /* } */
     494        // Test element only test
     495        bool active_element;
     496        this->element->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum);
     497        if(active_element){
     498                IssmDouble moulin_load,dt;
     499                this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
     500                parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     501                pe->values[0]=moulin_load*dt;
     502        }
     503        else pe->values[0]=0.0;
     504
     505
     506        //Test only input is sed
     507        /* pe->values[0]=0.0; */
     508
     509        //Clean up and return
    470510        return pe;
    471511}
Note: See TracChangeset for help on using the changeset viewer.