Changeset 23370


Ignore:
Timestamp:
10/02/18 10:00:01 (6 years ago)
Author:
bdef
Message:

BUG:partial fix, issues remain in Hydro

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

Legend:

Unmodified
Added
Removed

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

    r23368 r23370  
    607607                                //thickness[i] = old_thickness[i]*(1.0+opening-closing)/(1.0-opening+closing);
    608608
    609                                 /* if(element->nodes[i]->Sid()==2299){ */
    610                                 /*      printf("for node %i \n",element->nodes[i]->Sid()); */
    611                                 /*      printf(" old thickness is  %g \n",old_thickness[i]); */
    612                                 /*      printf(" new thickness is  %g \n",thickness[i]); */
    613                                 /*      printf(" closing is %g \n",((2.0*A*dt*pow(EPL_N,n))/(pow(n,n)))); */
    614                                 /*      printf(" opening is %g \n",((rho_water*gravity*epl_conductivity*EPLgrad2*dt)/(rho_ice*latentheat))); */
    615                                 /* } */
     609                                if(element->nodes[i]->Sid()==2299){
     610                                        printf("for node %i \n",element->nodes[i]->Sid());
     611                                        printf(" old thickness is  %g \n",old_thickness[i]);
     612                                        printf(" new thickness is  %g \n",thickness[i]);
     613                                }
    616614                                /*Take care of otherthikening*/
    617615                                if(thickness[i]>max_thick){
    … …  
    711709                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    712710                        /* If epl thickness gets under colapse thickness, close the layer */
    713                                 /* if(element->nodes[i]->Sid()==2299){ */
    714                                 /*      printf("thickness for maksing for node %i is %g \n",element->nodes[i]->Sid(),epl_thickness[i]);} */
    715711                        if(epl_thickness[i]<colapse_thick){
    716                                 /* if(element->nodes[i]->Sid()==2299){ */
    717                                 /*      printf("maksing node %i\n",element->nodes[i]->Sid());} */
     712                                if(element->nodes[i]->Sid()==2299){
     713                                        printf("masking node %i with thickness %f\n",element->nodes[i]->Sid(),epl_thickness[i]);}
    718714                                vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
    719715                                recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    720716                        }
    721                         /* //If epl head gets under base elevation, close the layer */
    722                         /* else if(eplhead[i]<(base[i]-1.0e-8)){ */
    723                         /*      vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); */
    724                         /*      recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); */
    725                         /* } */
    726717                }
    727718                /*If node is now closed bring its thickness back to initial*/

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

    r23366 r23370  
    2929
    3030        /*retrieve some parameters: */
     31        /* bool   issmb; */
     32        /* iomodel->FindConstant(&issmb,"md.transient.issmb"); */
    3133        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
    3234
    … …  
    6365                parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
    6466        }
     67        /* if(!issmb){ */
     68        /*      parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum)); */
     69        /* } */
     70
    6571  /*Requested outputs*/
    6672  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
    … …  
    7783        /*Fetch data needed: */
    7884        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
     85        bool   issmb;
     86        iomodel->FindConstant(&issmb,"md.transient.issmb");
    7987
    8088        /*Now, do we really want DC?*/
    … …  
    109117                iomodel->FetchDataToInput(elements,"md.initialization.epl_head",EplHeadHydrostepEnum);
    110118        }
     119        if(!issmb){
     120                iomodel->FetchDataToInput(elements,"md.smb.model",SmbEnum);
     121}
     122
    111123}/*}}}*/
    112124
    … …  
    638650        IssmDouble base_elev,prestep_head,water_sheet;
    639651        IssmDouble sediment_thickness       = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
    640         bool isthermal;
    641 
    642         element->FindParam(&isthermal,TransientIsthermalEnum);
     652
    643653        element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
    644654        SedTrans_input->GetInputValue(&FullLayer_transmitivity,gauss);

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

    r23335 r23370  
    188188
    189189        if(!issmb){
    190                 iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum);
     190                iomodel->FetchDataToInput(elements,"md.smb.model",SmbEnum);
    191191        }
    192192        if(stabilization==3){
    … …  
    629629}/*}}}*/
    630630void           MasstransportAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    631         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
     631        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    632632         * For node i, Bi can be expressed in the actual coordinate system
    633          * by: 
     633         * by:
    634634         *       Bi=[ N ]
    635635         *          [ N ]
    … …  
    657657}/*}}}*/
    658658void           MasstransportAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    659         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
     659        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    660660         * For node i, Bi' can be expressed in the actual coordinate system
    661          * by: 
     661         * by:
    662662         *       Bi_prime=[ dN/dx ]
    663663         *                [ dN/dy ]

  • issm/trunk-jpl/src/c/classes/Params/Parameters.cpp

    r23310 r23370  
    9292}
    9393/*}}}*/
    94 void Parameters::DeepEcho(void){/*{{{*/ 
     94void Parameters::DeepEcho(void){/*{{{*/
    9595        for(int i=0;i<NUMPARAMS;i++) {
    9696                if(this->params[i]) this->params[i]->DeepEcho();
    … …  
    9999}
    100100/*}}}*/
    101 void Parameters::Echo(void){/*{{{*/     
     101void Parameters::Echo(void){/*{{{*/
    102102        for(int i=0;i<NUMPARAMS;i++) {
    103103                if(this->params[i]) this->params[i]->Echo();
    … …  
    139139
    140140                        /*Recover enum of object first: */
    141                         MARSHALLING(obj_enum); 
     141                        MARSHALLING(obj_enum);
    142142
    143143                        if(obj_enum==DoubleParamEnum){
    … …  
    200200                                fileparam->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    201201                                delete fileparam;
    202                                 /* No need to add this object, the pointer is not valid             
     202                                /* No need to add this object, the pointer is not valid
    203203                                        The FemModel should reset all FileParams in the restart function */
    204204                        }
    … …  
    649649char* OptionsFromAnalysis(char** pouttoolkit,Parameters* parameters,int analysis_type){ /*{{{*/
    650650
    651         /* figure out ISSM options for current analysis, return a string. */ 
     651        /* figure out ISSM options for current analysis, return a string. */
    652652
    653653        /*output: */
    … …  
    667667
    668668        parameters->FindParam(&strings,&numanalyses,ToolkitsOptionsStringsEnum);
    669         parameters->FindParam(&toolkits,&dummy,ToolkitsTypesEnum); _assert_(dummy==numanalyses); 
    670         parameters->FindParam(&analyses,&dummy,ToolkitsOptionsAnalysesEnum); _assert_(dummy==numanalyses); 
     669        parameters->FindParam(&toolkits,&dummy,ToolkitsTypesEnum); _assert_(dummy==numanalyses);
     670        parameters->FindParam(&analyses,&dummy,ToolkitsOptionsAnalysesEnum); _assert_(dummy==numanalyses);
    671671
    672672        if(numanalyses==0)return NULL; //we did not find petsc options, don't bother.
    … …  
    710710        xDelete<int>(analyses);
    711711        return outstring;
    712 } 
     712}
    713713/*}}}*/
    714714void ToolkitsOptionsFromAnalysis(Parameters* parameters,int analysis_type){ /*{{{*/
    715715
    716716        /*!\file:  ToolkitsOptionsFromAnalysis.cpp
    717          * \brief: for each analysis, setup the issmoptions string. 
    718          * This is mainly for the case where we run our toolkits using petsc. In this case, we need to 
    719          * plug our toolkits options directly into the petsc options database. This is the case for each analysis type 
     717         * \brief: for each analysis, setup the issmoptions string.
     718         * This is mainly for the case where we run our toolkits using petsc. In this case, we need to
     719         * plug our toolkits options directly into the petsc options database. This is the case for each analysis type
    720720         * and parameters
    721          */ 
     721         */
    722722
    723723        char* options = NULL;
    … …  
    731731
    732732        #ifdef _HAVE_PETSC_
    733                 /*In case we are using PETSC, we do not rely on issmoptions. Instead, we dump issmoptions into the Petsc 
     733                /*In case we are using PETSC, we do not rely on issmoptions. Instead, we dump issmoptions into the Petsc
    734734                 * options database: */
    735735                #if (_PETSC_MINOR_>=7)

  • issm/trunk-jpl/src/c/cores/hydrology_core.cpp

    r23232 r23370  
    1414        /*Start profiler*/
    1515        femmodel->profiler->Start(HYDROLOGYCORE);
    16        
     16
    1717        /*intermediary*/
    1818        int          hydrology_model;
    … …  
    5050                /*intermediary: */
    5151                bool       isefficientlayer;
    52                 bool       isthermal;
    5352                int        step,hydroslices;
    5453                IssmDouble time,init_time,hydrotime,yts;
    … …  
    6261                femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
    6362
    64                 femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
    6563                /*first we exclude frozen nodes of the solved nodes*/
    6664                femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
    … …  
    153151                xDelete<char*>(requested_outputs);
    154152        }
    155        
     153
    156154        /*End profiler*/
    157155        femmodel->profiler->Stop(HYDROLOGYCORE);

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r23366 r23370  
    384384                }
    385385
    386                 /*duplicating smb position to have runoff value*/
     386                /*shifting smb position to have runoff value*/
    387387                if(issmb) smb_core(femmodel);
    388388
    … …  
    396396
    397397                /* from here on, prepare geometry for next time step*/
    398                 if(issmb) smb_core(femmodel);
     398                //if(issmb) smb_core(femmodel);
    399399
    400400                if(ismasstransport){
Note: See TracChangeset for help on using the changeset viewer.