source: issm/oecreview/Archive/23390-24306/ISSM-23543-23544.diff

Last change on this file was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 12.2 KB

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

     
    7474                                        int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum};
    7575                                        int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum};
    7676                                        int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
    77                                         femmodel->InitTransientOutputx(&stackedinput[0],4);
     77                                        //femmodel->InitTransientOutputx(&stackedinput[0],4);
     78                                        femmodel->InitMeanOutputx(&stackedinput[0],4);
    7879                                        while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
    7980                                                hydrostep+=1;
    8081                                                hydrotime+=hydrodt;
    … …  
    8788                                                /*Proceed now to heads computations*/
    8889                                                solutionsequence_hydro_nonlinear(femmodel);
    8990                                                /*If we have a sub-timestep we stack the variables here*/
    90                                                 femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
     91                                                //femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
     92                                                femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],4);
    9193                                        }
    92                                         femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
     94                                        //femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
     95                                        femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],4);
    9396                                }
    9497                                else{
    9598                                        int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum};
    9699                                        int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum};
    97100                                        int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum};
    98                                         femmodel->InitTransientOutputx(&stackedinput[0],2);
     101                                        //femmodel->InitTransientOutputx(&stackedinput[0],2);
     102                                        femmodel->InitMeanOutputx(&stackedinput[0],2);
    99103                                        while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
    100104                                                hydrotime+=hydrodt;
    101105                                                /*save preceding timestep*/
    … …  
    103107                                                /*Proceed now to heads computations*/
    104108                                                solutionsequence_hydro_nonlinear(femmodel);
    105109                                                /*If we have a sub-timestep we stack the variables here*/
    106                                                 femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
     110                                                //femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
     111                                                femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],2);
    107112                                        }
    108                                         femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
     113                                        //femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
     114                                        femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],2);
    109115                                }
    110116                        }
    111117                        else{

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

     
    48204820        if(VerboseSolution()) _printf0_("   Number of active nodes L2 Projection: "<< counter <<"\n");
    48214821}
    48224822/*}}}*/
    4823 void FemModel::InitTransientOutputx(int* input_enum,int numoutputs){ /*{{{*/
     4823void FemModel::InitTransientOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/
    48244824
    48254825  for(int i=0;i<numoutputs;i++){
    4826                 if(input_enum[i]<0){
     4826                if(stackedinput_enum[i]<0){
    48274827                        _error_("Can't deal with non enum fields for result Stack");
    48284828                }
    48294829                else{
    48304830                        for(int j=0;j<elements->Size();j++){
    4831                                 TransientInput* transient_input = new TransientInput(input_enum[i]);
     4831                                TransientInput* transient_input = new TransientInput(stackedinput_enum[i]);
    48324832                                /*Intermediaries*/
    48334833                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
    48344834                                //transient_input->Configure(element->parameters);
    … …  
    48494849                                /*Intermediaries*/
    48504850                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
    48514851                                TransientInput* stacking_input=NULL;
    4852                                 Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input);
    4853                                 Input* input_to_stack=element->GetInput(input_enum[i]); _assert_(input_to_stack);
     4852                                Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input); //this is the enum stack
    48544853                                stacking_input=dynamic_cast<TransientInput*>(input);
    48554854
    48564855                                int  numvertices = element->GetNumberOfVertices();
    48574856                                IssmDouble* N=xNew<IssmDouble>(numvertices);
    4858                                 element->GetInputListOnVertices(&N[0],input_enum[i]);
     4857                                element->GetInputListOnVertices(&N[0],input_enum[i]);   //this is the enum to stack
    48594858                                switch(element->ObjectEnum()){
    48604859                                case TriaEnum:
    48614860                                        stacking_input->AddTimeInput(new TriaInput(stackedinput_enum[i],&N[0],P1Enum),hydrotime);
    … …  
    48744873        }
    48754874}
    48764875/*}}}*/
    4877 void FemModel::AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/
     4876void FemModel::AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/
    48784877
    48794878  for(int i=0;i<numoutputs;i++){
    4880                 if(input_enum[i]<0){
     4879                if(stackedinput_enum[i]<0){
    48814880                        _error_("Can't deal with non enum fields for result Stack");
    48824881                }
    48834882                else{
    … …  
    48874886                                int  numvertices = element->GetNumberOfVertices();
    48884887                                IssmDouble* time_averaged=NULL;
    48894888
    4890                                 Input* input=element->inputs->GetInput(input_enum[i]); _assert_(input);
     4889                                Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input);
    48914890                                TransientInput* stacking_input=NULL;
    48924891                                stacking_input=dynamic_cast<TransientInput*>(input);
    48934892                                stacking_input->GetInputAverageOnTimes(&time_averaged,init_time);
    … …  
    48994898        }
    49004899}
    49014900/*}}}*/
     4901void FemModel::InitMeanOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/
     4902
     4903  for(int i=0;i<numoutputs;i++){
     4904                if(stackedinput_enum[i]<0){
     4905                        _error_("Can't deal with non enum fields for result Stack");
     4906                }
     4907                else{
     4908                        for(int j=0;j<elements->Size();j++){
     4909                                /*Intermediaries*/
     4910                                Element*   element            =xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     4911                                int        numvertices        =element->GetNumberOfVertices();
     4912                                IssmDouble zeros[numvertices] ={0.0};
     4913                                switch(element->ObjectEnum()){
     4914                                case TriaEnum:
     4915                                        element->inputs->AddInput(new TriaInput(stackedinput_enum[i],&zeros[0],P1Enum));
     4916                                        break;
     4917                                case PentaEnum:
     4918                                        element->inputs->AddInput(new PentaInput(stackedinput_enum[i],&zeros[0],P1Enum));
     4919                                        break;
     4920                                case TetraEnum:
     4921                                        element->inputs->AddInput(new TetraInput(stackedinput_enum[i],&zeros[0],P1Enum));
     4922                                        break;
     4923                                default: _error_("Not implemented yet");
     4924                                }
     4925                        }
     4926                }
     4927        }
     4928}
     4929/*}}}*/
     4930void FemModel::SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs){ /*{{{*/
     4931
     4932        //First get sub-timestep
     4933        IssmDouble hydrodt;
     4934        this->parameters->FindParam(&hydrodt,HydrologydtEnum);
     4935
     4936  for(int i=0;i<numoutputs;i++){
     4937                if(input_enum[i]<0){
     4938                        _error_("Can't deal with non enum fields for result Stack");
     4939                }
     4940                else{
     4941                        for(int j=0;j<elements->Size();j++){
     4942                                /*Intermediaries*/
     4943                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     4944                                int  numvertices = element->GetNumberOfVertices();
     4945                                IssmDouble* values_to_add=xNew<IssmDouble>(numvertices);
     4946                                IssmDouble* existing_values=xNew<IssmDouble>(numvertices);
     4947                                element->GetInputListOnVertices(&values_to_add[0],input_enum[i]); //those are the values to add
     4948                                element->GetInputListOnVertices(&existing_values[0],stackedinput_enum[i]); //those are the values to add
     4949                                for(int k=0;k<numvertices;k++){
     4950                                        existing_values[k]+=values_to_add[k]*hydrodt;
     4951                                }
     4952                                element->AddInput(stackedinput_enum[i],&existing_values[0],P1Enum);
     4953                                xDelete<IssmDouble>(existing_values);
     4954                                xDelete<IssmDouble>(values_to_add);
     4955                        }
     4956                }
     4957        }
     4958}
     4959/*}}}*/
     4960void FemModel::AverageSumOutputx(int* stackedinput_enum,int* averagedinput_enum,int numoutputs){ /*{{{*/
     4961
     4962        //First get timestep
     4963        IssmDouble maindt;
     4964        this->parameters->FindParam(&maindt,TimesteppingTimeStepEnum);
     4965  for(int i=0;i<numoutputs;i++){
     4966                if(stackedinput_enum[i]<0){
     4967                        _error_("Can't deal with non enum fields for result Stack");
     4968                }
     4969                else{
     4970                        for(int j=0;j<elements->Size();j++){
     4971                                /*Intermediaries*/
     4972                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     4973                                int  numvertices = element->GetNumberOfVertices();
     4974                                IssmDouble* time_averaged=xNew<IssmDouble>(numvertices);
     4975                                IssmDouble* existing_values=xNew<IssmDouble>(numvertices);
     4976                                element->GetInputListOnVertices(&existing_values[0],stackedinput_enum[i]); //those are the values to add
     4977
     4978                                for(int k=0;k<numvertices;k++){
     4979                                        time_averaged[k]=existing_values[k]/maindt;
     4980                                }
     4981
     4982                                element->AddInput(averagedinput_enum[i],&time_averaged[0],P1Enum);
     4983                                xDelete<IssmDouble>(time_averaged);
     4984                                xDelete<IssmDouble>(existing_values);
     4985                        }
     4986                }
     4987        }
     4988}
     4989/*}}}*/
    49024990#ifdef _HAVE_JAVASCRIPT_
    49034991FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
    49044992        /*configuration: */

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

     
    140140                void Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y);
    141141                #endif
    142142                #ifdef _HAVE_ESA_
    143                 void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy); 
    144                 void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); 
     143                void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy);
     144                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    145145                #endif
    146146                #ifdef _HAVE_SEALEVELRISE_
    147147                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
    148148                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
    149149                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    150                 void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); 
     150                void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
    151151                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
    152152                #endif
    153153                void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);
    … …  
    158158                void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
    159159                void InitTransientOutputx(int* input_enum, int numoutputs);
    160160                void StackTransientOutputx(int* input_enum,int* stackedinput_enum, IssmDouble hydrotime, int numoutputs);
    161                 void AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble hydrotime,int numoutputs);
     161                void AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs);
     162                void InitMeanOutputx(int* stackedinput_enum,int numoutputs);
     163                void SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs);
     164                void AverageSumOutputx(int* stackedinput_enum,int* averagedinput_enum,int numoutputs);
    162165                void UpdateConstraintsx(void);
    163166                int  UpdateVertexPositionsx(void);
    164167
Note: See TracBrowser for help on using the repository browser.