Changeset 22277


Ignore:
Timestamp:
Nov 28, 2017, 9:35:29 AM (9 years ago)
Author:
bdef
Message:

NEW: Implementing substep framework in hydrologydc

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

Legend:

Unmodified
Added
Removed

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

    r22236 r22277  
    223223        /*Retrieve all inputs and parameters*/
    224224        basalelement->GetVerticesCoordinates(&xyz_list);
    225         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     225        //basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     226        basalelement ->FindParam(&dt,HydrologydcHydrodtEnum);
    226227
    227228        Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input);
    … …  
    332333        /*Retrieve all inputs and parameters*/
    333334        basalelement->GetVerticesCoordinates(&xyz_list);
    334         basalelement->FindParam(&dt,TimesteppingTimeStepEnum); 
     335        //basalelement->FindParam(&dt,TimesteppingTimeStepEnum);       
     336        basalelement ->FindParam(&dt,HydrologydcHydrodtEnum);
    335337
    336338        Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input);
    … …  
    555557                Input*  active_element_input=element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);               
    556558                active_element_input->GetInputValue(&active_element);
    557                 element->FindParam(&dt,TimesteppingTimeStepEnum);
     559                //element->FindParam(&dt,TimesteppingTimeStepEnum);
     560                element ->FindParam(&dt,HydrologydcHydrodtEnum);
    558561       
    559562                /*For now, assuming just one way to compute EPL thickness*/
    … …  
    725728        basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType());
    726729
    727         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     730        if(domaintype!=Domain2DhorizontalEnum){
     731                basalelement->DeleteMaterials();
     732                delete basalelement;
     733        }
    728734        xDelete<IssmDouble>(epl_thickness);
    729735        xDelete<IssmDouble>(old_active);

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

    r21964 r22277  
    1919        int         penalty_lock;
    2020        int         hydro_maxiter;
     21        int         hydroslices;
     22        int         numoutputs;
    2123        bool        isefficientlayer;
    2224        IssmDouble  penalty_factor;
    … …  
    2426        IssmDouble  leakagefactor;
    2527        IssmDouble  sedimentlimit;
     28        char**      requestedoutputs = NULL;
    2629
    2730        /*retrieve some parameters: */
    … …  
    3336        iomodel->FetchData(&sedimentlimit_flag, "md.hydrology.sedimentlimit_flag" );
    3437        iomodel->FetchData(&transfer_flag,      "md.hydrology.transfer_flag" );
    35         iomodel->FetchData(&unconfined_flag,      "md.hydrology.unconfined_flag" );
     38        iomodel->FetchData(&unconfined_flag,    "md.hydrology.unconfined_flag" );
     39        iomodel->FetchData(&penalty_lock,       "md.hydrology.penalty_lock" );
     40        iomodel->FetchData(&hydro_maxiter,      "md.hydrology.max_iter" );
     41        iomodel->FetchData(&hydroslices,        "md.hydrology.steps_per_step");
     42        iomodel->FetchData(&isefficientlayer,   "md.hydrology.isefficientlayer");
    3643        iomodel->FetchData(&penalty_factor,     "md.hydrology.penalty_factor" );
    37         iomodel->FetchData(&isefficientlayer,   "md.hydrology.isefficientlayer");
    38         iomodel->FetchData(&hydro_maxiter,      "md.hydrology.max_iter" );
    39         iomodel->FetchData(&penalty_lock,       "md.hydrology.penalty_lock" );
    4044        iomodel->FetchData(&rel_tol,            "md.hydrology.rel_tol" );
    4145
    … …  
    4650        parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
    4751        parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
     52        parameters->AddObject(new IntParam(HydrologydcStepsPerStepEnum,hydroslices));
     53
    4854        parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
    4955        parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
    … …  
    5763                parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
    5864        }
     65  /*Requested outputs*/
     66  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
     67  parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
     68  if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
     69  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
    5970}/*}}}*/
    6071
    … …  
    96107        if(isefficientlayer){
    97108                iomodel->FetchDataToInput(elements,"md.hydrology.mask_eplactive_node",HydrologydcMaskEplactiveNodeEnum);
     109                iomodel->FetchDataToInput(elements,"md.initialization.epl_head",EplHeadEnum);
    98110        }
    99111}/*}}}*/
    … …  
    209221        /*Retrieve all inputs and parameters*/
    210222        basalelement ->GetVerticesCoordinates(&xyz_list);
    211         basalelement ->FindParam(&dt,TimesteppingTimeStepEnum);
     223        //basalelement ->FindParam(&dt,TimesteppingTimeStepEnum);
     224        basalelement ->FindParam(&dt,HydrologydcHydrodtEnum);
    212225        basalelement ->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    213226        Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
    214227        Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum);
    215228        Input* base_input     = basalelement->GetInput(BaseEnum);
     229        Input* old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);                  _assert_(old_wh_input);
    216230
    217231        /*Transfer related Inputs*/
    … …  
    232246                /*Diffusivity*/
    233247                D_scalar=sediment_transmitivity*gauss->weight*Jdet;
     248                //D_scalar=gauss->weight*Jdet;
    234249                if(dt!=0.) D_scalar=D_scalar*dt;
    235250                D[0][0]=D_scalar;
    … …  
    245260                        basalelement->NodalFunctions(&basis[0],gauss);
    246261                        D_scalar=sediment_storing*gauss->weight*Jdet;
     262                        //D_scalar=(sediment_storing/sediment_transmitivity)*gauss->weight*Jdet;
    247263                        TripleMultiply(basis,numnodes,1,0,
    248264                                                                                 &D_scalar,1,1,0,
    … …  
    257273                                        basalelement->NodalFunctions(&basis[0],gauss);
    258274                                        D_scalar=dt*transfer*gauss->weight*Jdet;
     275                                        //D_scalar=dt*(transfer/sediment_transmitivity)*gauss->weight*Jdet;
    259276                                        TripleMultiply(basis,numnodes,1,0,
    260277                                                                                                 &D_scalar,1,1,0,
    … …  
    270287        xDelete<IssmDouble>(basis);
    271288        delete gauss;
    272         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     289        if(domaintype!=Domain2DhorizontalEnum){
     290                basalelement->DeleteMaterials();
     291                delete basalelement;
     292        }
    273293        return Ke;
    274294}/*}}}*/
    … …  
    296316        bool       active_element,isefficientlayer;
    297317        IssmDouble dt,scalar,sediment_storing;
    298         IssmDouble water_head;
     318        IssmDouble water_head,sediment_transmitivity;
    299319        IssmDouble water_load,transfer;
    300320        IssmDouble Jdet;
    … …  
    313333        /*Retrieve all inputs and parameters*/
    314334        basalelement->GetVerticesCoordinates(&xyz_list);
    315         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     335        //basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     336        basalelement ->FindParam(&dt,HydrologydcHydrodtEnum);
    316337        basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    317338
    … …  
    320341        Input* base_input                 = basalelement->GetInput(BaseEnum);
    321342        Input* water_input        = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
     343        Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
    322344
    323345        if(dt!= 0.){
    … …  
    335357                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    336358                basalelement->NodalFunctions(basis,gauss);
     359                sediment_transmitivity = SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
    337360
    338361                /*Loading term*/
    … …  
    340363                        water_input->GetInputValue(&water_load,gauss);
    341364                        scalar = Jdet*gauss->weight*(water_load);
     365                        //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
    342366                        if(dt!=0.) scalar = scalar*dt;
    343367                        for(int i=0;i<numnodes;i++){
    … …  
    351375                                water_input->GetInputValue(&water_load,gauss);
    352376                                scalar = Jdet*gauss->weight*(water_load);
     377                                //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
    353378                                if(dt!=0.) scalar = scalar*dt;
    354379                                for(int i=0;i<numnodes;i++){
    … …  
    372397                                }
    373398                                scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer));
     399                                //scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer))/sediment_transmitivity;
    374400                                for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    375401                        }
    376402                        else{
    377403                                scalar = Jdet*gauss->weight*(water_head*sediment_storing);
     404                                //scalar = Jdet*gauss->weight*(water_head*sediment_storing)/sediment_transmitivity;
    378405                                for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    379406                        }
    … …  
    384411        xDelete<IssmDouble>(basis);
    385412        delete gauss;
    386         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     413        if(domaintype!=Domain2DhorizontalEnum){
     414                basalelement->DeleteMaterials();
     415                delete basalelement;
     416        }
    387417        return pe;
    388418}/*}}}*/
    … …  
    470500                IssmDouble* thickness = xNew<IssmDouble>(numnodes);
    471501                IssmDouble* base      = xNew<IssmDouble>(numnodes);
     502                IssmDouble* transmitivity = xNew<IssmDouble>(numnodes);
    472503
    473504                basalelement->FindParam(&kmax,HydrologySedimentKmaxEnum);
    … …  
    479510                basalelement->GetInputListOnVertices(thickness,ThicknessEnum);
    480511                basalelement->GetInputListOnVertices(base,BaseEnum);
    481 
     512                basalelement->GetInputListOnVertices(transmitivity,HydrologydcSedimentTransmitivityEnum);
    482513                kappa=kmax*pow(10.,penalty_factor);
    483514
    … …  
    486517                        if(values[i]>h_max) {
    487518                                residual[i] = kappa*(values[i]-h_max);
     519                                //residual[i] = kappa*(values[i]-h_max)*transmitivity[i];
    488520                        }
    489521                        else{
    … …  
    495527                }
    496528                xDelete<IssmDouble>(thickness);
     529                xDelete<IssmDouble>(transmitivity);
    497530                xDelete<IssmDouble>(base);
    498531        }
    … …  
    500533        /*Add input to the element: */
    501534        element->AddBasalInput(SedimentHeadEnum,values,P1Enum);
     535        element->AddBasalInput(EffectivePressureEnum,pressure,P1Enum);
    502536        element->AddBasalInput(SedimentHeadResidualEnum,residual,P1Enum);
    503         element->AddBasalInput(EffectivePressureEnum,pressure,P1Enum);
    504537
    505538        /*Free ressources:*/
    … …  
    508541        xDelete<IssmDouble>(pressure);
    509542        xDelete<int>(doflist);
    510         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     543        if(domaintype!=Domain2DhorizontalEnum){
     544                basalelement->DeleteMaterials();
     545                delete basalelement;
     546        }
    511547}/*}}}*/
    512548
    … …  
    520556        IssmDouble expfac;
    521557        IssmDouble sediment_storing;
    522         IssmDouble storing;
     558        IssmDouble storing,yield;
    523559        IssmDouble base_elev,prestep_head,water_sheet;
    524560        IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    … …  
    538574                water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
    539575
     576                /* if (water_sheet<sediment_thickness){ */
     577                /*      sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); */
     578                /* } */
     579                /* else{ */
     580                /*      sediment_storing=sediment_porosity; */
     581                /* } */
    540582                storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
    541        
    542                 //      using logistic function for heavyside approximation
     583                //using logistic function for heavyside approximation
    543584                expfac=10.;
    544                 sediment_storing=sediment_porosity+(storing-sediment_porosity)/(1+exp(-2*expfac*(water_sheet-0.99*sediment_thickness)));
     585                yield=sediment_porosity;
     586                sediment_storing=yield+(storing-yield)/(1+exp(-2*expfac*(water_sheet-0.99*sediment_thickness)));
    545587                break;
    546588        default:
    … …  
    568610                sed_head_input->GetInputValue(&prestep_head,gauss);
    569611                water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
    570 
    571                 if (water_sheet<=sediment_thickness){
    572                         sediment_transmitivity=FullLayer_transmitivity*water_sheet/sediment_thickness;
    573                 }
    574                 else{
    575                         sediment_transmitivity=FullLayer_transmitivity;
    576                 }
     612               
     613                //min definition of the if test
     614                sediment_transmitivity=FullLayer_transmitivity*min(water_sheet/sediment_thickness,1.);
     615                if (sediment_transmitivity==0){
     616                        sediment_transmitivity=1.0e-20;
     617                }
     618
    577619                break;
    578620        default:
    … …  
    632674        case 1:
    633675                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    634                 transfer=leakage;
     676                transfer=+leakage;
    635677                break;
    636678        default:
    … …  
    656698                _assert_(epl_head_input);
    657699                epl_head_input->GetInputValue(&epl_head,gauss);
     700                if (element->Id()==42){
     701                        printf("epl head in sed Pvec transfer is %f\n",epl_head);
     702                }
    658703                element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    659                 transfer=epl_head*leakage;
     704                transfer=+epl_head*leakage;
    660705                break;
    661706        default:

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r22266 r22277  
    17731773                                name==HydrologydcMaskEplactiveNodeEnum ||
    17741774                                name==HydrologyHeadEnum ||
    1775                  name==HydrologyHeadOldEnum ||         
     1775                                name==HydrologyHeadOldEnum ||           
    17761776                                name==StressbalanceConvergenceNumStepsEnum ||
    17771777                                name==MeshVertexonbaseEnum ||

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r22276 r22277  
    44584458}
    44594459/*}}}*/
     4460void FemModel::InitTransientOutputx(int* input_enum,int numoutputs){ /*{{{*/
     4461
     4462  for(int i=0;i<numoutputs;i++){
     4463                if(input_enum[i]<0){
     4464                        _error_("Can't deal with non enum fields for result Stack");
     4465                }
     4466                else{
     4467                        for(int j=0;j<elements->Size();j++){
     4468                                TransientInput* transient_input = new TransientInput(input_enum[i]);
     4469                                /*Intermediaries*/
     4470                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     4471                                transient_input->Configure(element->parameters);
     4472                                element->inputs->AddInput(transient_input);
     4473                        }
     4474                }
     4475        }
     4476}
     4477/*}}}*/
     4478void FemModel::StackTransientOutputx(int* input_enum,IssmDouble hydrotime,int numoutputs){ /*{{{*/
     4479
     4480  for(int i=0;i<numoutputs;i++){
     4481                if(input_enum[i]<0){
     4482                        _error_("Can't deal with non enum fields for result Stack");
     4483                }
     4484                else{
     4485                        for(int j=0;j<elements->Size();j++){
     4486                                /*Intermediaries*/
     4487                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     4488                                TransientInput* stacking_input=NULL;
     4489                                Input* input=element->inputs->GetInput(HydrologydcStackedNEnum); _assert_(input);
     4490                                Input* input_to_stack=element->GetInput(input_enum[i]); _assert_(input_to_stack);
     4491                                stacking_input=dynamic_cast<TransientInput*>(input);
     4492
     4493                                int  numvertices = element->GetNumberOfVertices();
     4494                                IssmDouble* N=xNew<IssmDouble>(numvertices);
     4495                                element->GetInputListOnVertices(&N[0],input_enum[i]);
     4496                                switch(element->ObjectEnum()){
     4497                                case TriaEnum:
     4498                                        stacking_input->AddTimeInput(new TriaInput(HydrologydcStackedNEnum,&N[0],P1Enum),hydrotime);
     4499                                        break;
     4500                                case PentaEnum:
     4501                                        stacking_input->AddTimeInput(new PentaInput(HydrologydcStackedNEnum,&N[0],P1Enum),hydrotime);
     4502                                        break;
     4503                                case TetraEnum:
     4504                                        stacking_input->AddTimeInput(new TetraInput(HydrologydcStackedNEnum,&N[0],P1Enum),hydrotime);
     4505                                        break;
     4506                                default: _error_("Not implemented yet");
     4507                                }
     4508                                xDelete<IssmDouble>(N);
     4509                        }
     4510                }
     4511        }
     4512}
     4513/*}}}*/
     4514void FemModel::AverageTransientOutputx(int* input_enum,IssmDouble init_time,int numoutputs){ /*{{{*/
     4515        IssmDouble* time_averaged=NULL;
     4516  for(int i=0;i<numoutputs;i++){
     4517                if(input_enum[i]<0){
     4518                        _error_("Can't deal with non enum fields for result Stack");
     4519                }
     4520                else{
     4521                        for(int j=0;j<elements->Size();j++){
     4522                                /*Intermediaries*/
     4523                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     4524                                int  numvertices = element->GetNumberOfVertices();
     4525                                //xNew<IssmDouble>(numvertices);
     4526
     4527                                Input* input=element->inputs->GetInput(input_enum[i]); _assert_(input);
     4528                                TransientInput* stacking_input=NULL;
     4529                                stacking_input=dynamic_cast<TransientInput*>(input);
     4530                                stacking_input->GetInputAverageOnTimes(&time_averaged, init_time);
     4531
     4532                                element->AddInput(TimeAverageEffectivePressureEnum,time_averaged,P1Enum);
     4533                        }
     4534                }
     4535        }
     4536}
     4537/*}}}*/
    44604538#ifdef _HAVE_JAVASCRIPT_
    44614539FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r22266 r22277  
    149149                void UpdateConstraintsExtrudeFromTopx();
    150150                void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
     151                void InitTransientOutputx(int* input_enum, int numoutputs);
     152                void StackTransientOutputx(int* input_enum, IssmDouble hydrotime, int numoutputs);
     153                void AverageTransientOutputx(int* input_enum,IssmDouble hydrotime,int numoutputs);
    151154                void UpdateConstraintsx(void);
    152155                int  UpdateVertexPositionsx(void);

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp

    r21920 r22277  
    272272
    273273        int         i;
    274         IssmDouble *times  = NULL;
    275         IssmDouble *values = NULL;
     274        IssmDouble* times  = NULL;
     275        IssmDouble* values = NULL;
    276276        int         numsteps;
    277277        bool        iscurrenttime_included = false;
    … …  
    291291
    292292        for(i=0;i<numsteps;i++){
    293 
    294293                if((iscurrenttime_included==false) && (i==(numsteps-1))){
    295 
    296294                        /*Retrieve interpolated values for current time step: */
    297295                        Input* input=GetTimeInput(currenttime);
    … …  
    311309}
    312310/*}}}*/
     311void TransientInput::GetInputAverageOnTimes(IssmDouble** pvalues, IssmDouble init_time){/*{{{*/
     312
     313        IssmDouble  time;
     314        IssmDouble  preceding_time;
     315        IssmDouble  values[3] = {0.0,0.0,0.0};
     316
     317        for(int i=0;i<this->numtimesteps;i++){
     318                time                                             = this->timesteps[i];
     319                if(i==0){
     320                        preceding_time = init_time;
     321                }
     322                else{
     323                        preceding_time = this->timesteps[i-1];
     324                }
     325                TriaInput*      input            = (TriaInput*)this->inputs->GetObjectByOffset(i); _assert_(input);
     326                int                                     numnodes = input->NumberofNodes(input->interpolation_type);
     327                for(int j=0;j<numnodes;j++){
     328                        values[j]+=(input->values[j]*(time-preceding_time));
     329                }
     330        }
     331        for(int j       =       0;j<3;j++){
     332                values[j]/=(this->timesteps[this->numtimesteps-1]-init_time);
     333                //printf("final value is %e\n",values[j]);
     334        }
     335        *pvalues=values;
     336}
     337/*}}}*/
    313338
    314339/*Intermediary*/

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h

    r21974 r22277  
    6464                void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
    6565                void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime);
     66                void GetInputAverageOnTimes(IssmDouble** pvalues, IssmDouble init_time);
    6667                int  GetInputInterpolationType(){_error_("not implemented yet!");};
    6768                IssmDouble  GetTimeByOffset(int offset);

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp

    r21974 r22277  
    223223}
    224224/*}}}*/
     225
    225226void TriaInput::GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){/*{{{*/
    226227

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

    r22185 r22277  
    1414        /*intermediary*/
    1515        int  hydrology_model;
     16        int  solution_type;
    1617        bool save_results;
    1718        bool modify_loads=true;
    … …  
    2122        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    2223        femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
    23        
     24        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);       
     25
    2426        if(VerboseSolution()) _printf0_("   computing water heads\n");
    2527                       
    … …  
    5254        /*Using the double continuum model*/
    5355        else if (hydrology_model==HydrologydcEnum){
    54                 InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum);
    55                 femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    56                 if (isefficientlayer){
    57                         InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
    58                         InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum);
     56                /*intermediary: */
     57                int        step,hydroslices;
     58                int        numoutputs;
     59                char       **requested_outputs = NULL;
     60                IssmDouble time,init_time,hydrotime,yts;
     61                IssmDouble dt,hydrodt;
     62
     63                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     64                femmodel->parameters->FindParam(&step,StepEnum);
     65                femmodel->parameters->FindParam(&time,TimeEnum);
     66                femmodel->parameters->FindParam(&hydroslices,HydrologydcStepsPerStepEnum);
     67                femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
     68                init_time = time-dt; //getting the time back to the start of the timestep
     69                hydrotime=init_time;
     70                hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider
     71                femmodel->parameters->AddObject(new DoubleParam(HydrologydcHydrodtEnum,hydrodt));
     72                femmodel->parameters->FindParam(&numoutputs,HydrologyNumRequestedOutputsEnum);
     73                if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,HydrologyRequestedOutputsEnum);
     74
     75                if(dt>0){
     76                        if(hydroslices>1){
     77                                int trans_input[1]={HydrologydcStackedNEnum};
     78                                femmodel->InitTransientOutputx(&trans_input[0],1);
     79                        }
     80                        while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
     81                                hydrotime+=hydrodt;
     82                                InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum);
     83                                femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     84                                if (isefficientlayer){
     85                                        InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
     86                                        InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum);
     87                                }
     88                                /*Proceed now to heads computations*/
     89                                solutionsequence_hydro_nonlinear(femmodel);
     90                                if (hydroslices>1){
     91                                        int output[1]={EffectivePressureEnum};
     92                                        femmodel->StackTransientOutputx(&output[0],hydrotime,1);
     93                                }
     94                        }
     95                        if(hydroslices>1){
     96                                int output[1]={HydrologydcStackedNEnum};
     97                                femmodel->AverageTransientOutputx(&output[0],init_time,1);
     98                        }               
    5999                }
    60                
    61                 /*Proceed now to heads computations*/
    62                 solutionsequence_hydro_nonlinear(femmodel);
    63 
     100                else{
     101                        InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum);
     102                        femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     103                        if (isefficientlayer){
     104                                InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
     105                                InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum);
     106                        }
     107                        /*Proceed now to heads computations*/
     108                        solutionsequence_hydro_nonlinear(femmodel);
     109                }
    64110                if(save_results){
    65111                        if(VerboseSolution()) _printf0_("   saving results \n");
    66                         if(isefficientlayer){
    67                                 int outputs[9] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,HydrologydcMaskEplactiveEltEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum,EffectivePressureEnum};
    68                                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],9);
    69                         }
    70                         else{
    71                                 int outputs[3] = {SedimentHeadEnum,SedimentHeadResidualEnum,EffectivePressureEnum};
    72                                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
    73                         }
    74                         /*unload results*/
    75                         if(VerboseSolution()) _printf0_("   saving temporary results\n");
    76                         OutputResultsx(femmodel);
     112                        femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
     113                        /*unload results removed 23/11 needs checking*/
     114                        /* if(VerboseSolution()) _printf0_("   saving temporary results\n"); */
     115                        /* OutputResultsx(femmodel); */
     116                }
     117                /*Free ressources:*/   
     118                if(numoutputs){
     119                        for(int i=0;i<numoutputs;i++){
     120                                xDelete<char>(requested_outputs[i]);
     121                        }
     122                        xDelete<char*>(requested_outputs);
    77123                }
    78124        }
     125       
    79126
    80127        else if (hydrology_model==HydrologysommersEnum){       

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r22276 r22277  
    131131        HydrologyshreveStabilizationEnum,
    132132        HydrologydcEnum,
     133        HydrologyNumRequestedOutputsEnum,
     134        HydrologyRequestedOutputsEnum,
     135        HydrologydcHydrodtEnum,
     136        HydrologydcStepsPerStepEnum,
    133137        SedimentHeadEnum,
    134138        SedimentHeadOldEnum,
    135139        SedimentHeadResidualEnum,
    136140        EffectivePressureEnum,
     141        TimeAverageEffectivePressureEnum,
    137142        EplHeadEnum,
    138143        EplHeadOldEnum,
    … …  
    140145        EplHeadSlopeYEnum,
    141146        EplZigZagCounterEnum,
     147        MeanEffectivePressureEnum,
     148        HydrologydcStackedNEnum,
    142149        HydrologydcMaxIterEnum,
    143150        HydrologydcRelTolEnum,
    … …  
    182189        HydrologyReynoldsEnum,
    183190        HydrologyNeumannfluxEnum,
    184    HydrologyRelaxationEnum,
     191        HydrologyRelaxationEnum,
    185192        HydrologyBasalFluxEnum,
    186193        HydrologyStorageEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r22276 r22277  
    137137                case HydrologyshreveStabilizationEnum : return "HydrologyshreveStabilization";
    138138                case HydrologydcEnum : return "Hydrologydc";
     139                case HydrologyNumRequestedOutputsEnum : return "HydrologyNumRequestedOutputs";
     140                case HydrologyRequestedOutputsEnum : return "HydrologyRequestedOutputs";
     141                case HydrologydcHydrodtEnum : return "HydrologydcHydrodt";
     142                case HydrologydcStepsPerStepEnum : return "HydrologydcStepsPerStep";
    139143                case SedimentHeadEnum : return "SedimentHead";
    140144                case SedimentHeadOldEnum : return "SedimentHeadOld";
    141145                case SedimentHeadResidualEnum : return "SedimentHeadResidual";
    142146                case EffectivePressureEnum : return "EffectivePressure";
     147                case TimeAverageEffectivePressureEnum : return "TimeAverageEffectivePressure";
    143148                case EplHeadEnum : return "EplHead";
    144149                case EplHeadOldEnum : return "EplHeadOld";
    … …  
    146151                case EplHeadSlopeYEnum : return "EplHeadSlopeY";
    147152                case EplZigZagCounterEnum : return "EplZigZagCounter";
     153                case MeanEffectivePressureEnum : return "MeanEffectivePressure";
     154                case HydrologydcStackedNEnum : return "HydrologydcStackedN";
    148155                case HydrologydcMaxIterEnum : return "HydrologydcMaxIter";
    149156                case HydrologydcRelTolEnum : return "HydrologydcRelTol";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r22276 r22277  
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
     142              if (strcmp(name,"HydrologyNumRequestedOutputs")==0) return HydrologyNumRequestedOutputsEnum;
     143              else if (strcmp(name,"HydrologyRequestedOutputs")==0) return HydrologyRequestedOutputsEnum;
     144              else if (strcmp(name,"HydrologydcHydrodt")==0) return HydrologydcHydrodtEnum;
     145              else if (strcmp(name,"HydrologydcStepsPerStep")==0) return HydrologydcStepsPerStepEnum;
     146              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    143147              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    144148              else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum;
    145149              else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
     150              else if (strcmp(name,"TimeAverageEffectivePressure")==0) return TimeAverageEffectivePressureEnum;
    146151              else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
    147152              else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
    … …  
    149154              else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
    150155              else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
     156              else if (strcmp(name,"MeanEffectivePressure")==0) return MeanEffectivePressureEnum;
     157              else if (strcmp(name,"HydrologydcStackedN")==0) return HydrologydcStackedNEnum;
    151158              else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum;
    152159              else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
    … …  
    253260              else if (strcmp(name,"DamageHealing")==0) return DamageHealingEnum;
    254261              else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum;
    255               else if (strcmp(name,"DamageKappa")==0) return DamageKappaEnum;
     262         else stage=3;
     263   }
     264   if(stage==3){
     265              if (strcmp(name,"DamageKappa")==0) return DamageKappaEnum;
    256266              else if (strcmp(name,"DamageStabilization")==0) return DamageStabilizationEnum;
    257267              else if (strcmp(name,"DamageMaxiter")==0) return DamageMaxiterEnum;
    … …  
    260270              else if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum;
    261271              else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum;
    262          else stage=3;
    263    }
    264    if(stage==3){
    265               if (strcmp(name,"Damage")==0) return DamageEnum;
     272              else if (strcmp(name,"Damage")==0) return DamageEnum;
    266273              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    267274              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
    … …  
    376383              else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
    377384              else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
    378               else if (strcmp(name,"TimesteppingInterpForcings")==0) return TimesteppingInterpForcingsEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"TimesteppingInterpForcings")==0) return TimesteppingInterpForcingsEnum;
    379389              else if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum;
    380390              else if (strcmp(name,"TransientIscoupler")==0) return TransientIscouplerEnum;
    … …  
    383393              else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
    384394              else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
     395              else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
    389396              else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum;
    390397              else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
    … …  
    499506              else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
    500507              else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
    501               else if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum;
    502512              else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
    503513              else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
    … …  
    506516              else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
    507517              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
     518              else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
    512519              else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
    513520              else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
    … …  
    622629              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
    623630              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
    624               else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
    625635              else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
    626636              else if (strcmp(name,"Outputdefinition1")==0) return Outputdefinition1Enum;
    … …  
    629639              else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
    630640              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
     641              else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
    635642              else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
    636643              else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
    … …  
    745752              else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
    746753              else if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum;
    747               else if (strcmp(name,"SubelementMigration3")==0) return SubelementMigration3Enum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"SubelementMigration3")==0) return SubelementMigration3Enum;
    748758              else if (strcmp(name,"Contact")==0) return ContactEnum;
    749759              else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
    … …  
    752762              else if (strcmp(name,"Colinear")==0) return ColinearEnum;
    753763              else if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"Fset")==0) return FsetEnum;
     764              else if (strcmp(name,"Fset")==0) return FsetEnum;
    758765              else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
    759766              else if (strcmp(name,"Gradient2")==0) return Gradient2Enum;
    … …  
    868875              else if (strcmp(name,"AmrRestart")==0) return AmrRestartEnum;
    869876              else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
    870               else if (strcmp(name,"AmrLevelMax")==0) return AmrLevelMaxEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"AmrLevelMax")==0) return AmrLevelMaxEnum;
    871881              else if (strcmp(name,"AmrLag")==0) return AmrLagEnum;
    872882              else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
    … …  
    875885              else if (strcmp(name,"AmrField")==0) return AmrFieldEnum;
    876886              else if (strcmp(name,"AmrErr")==0) return AmrErrEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"AmrKeepMetric")==0) return AmrKeepMetricEnum;
     887              else if (strcmp(name,"AmrKeepMetric")==0) return AmrKeepMetricEnum;
    881888              else if (strcmp(name,"AmrGradation")==0) return AmrGradationEnum;
    882889              else if (strcmp(name,"AmrGroundingLineResolution")==0) return AmrGroundingLineResolutionEnum;
    … …  
    991998              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
    992999              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    993               else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    9941004              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    9951005              else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
    … …  
    9981008              else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    9991009              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum;
     1010              else if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum;
    10041011              else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    10051012              else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
Note: See TracChangeset for help on using the changeset viewer.