Changeset 24565


Ignore:
Timestamp:
02/14/20 05:45:36 (5 years ago)
Author:
bdef
Message:

BUG:fix for the coupling on different time steps

Location:
issm/trunk-jpl
Files:
20 edited

Legend:

Unmodified
Added
Removed

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

    r24529 r24565  
    335335        int        smb_model;
    336336        int        smbsubstepping;
     337        int        hydrologysubstepping;
    337338        IssmDouble dt,scalar,sediment_storing;
    338339        IssmDouble water_head,sediment_transmitivity;
    … …  
    340341        IssmDouble Jdet;
    341342
    342         IssmDouble      *xyz_list             = NULL;
    343         Input2          *active_element_input = NULL;
    344         Input2          *old_wh_input         = NULL;
    345         Input2          *dummy_input          = NULL;
    346         TransientInput2 *surface_runoff_input = NULL;
     343        IssmDouble *xyz_list             = NULL;
     344        Input2     *active_element_input = NULL;
     345        Input2     *old_wh_input         = NULL;
     346        Input2     *dummy_input          = NULL;
     347        Input2    *surface_runoff_input = NULL;
    347348
    348349        /*Fetch number of nodes and dof for this finite element*/
    … …  
    370371
    371372        if(dt!= 0.){
    372                 old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum);                  _assert_(old_wh_input);
     373                old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input);
    373374        }
    374375        if(smb_model==SMBgradientscomponentsEnum){
    375376                basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum);
    376                 if(smbsubstepping>1) {
    377                         dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum); _assert_(dummy_input);
    378                 }
    379                 else {
     377                basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum);
     378
     379                if(smbsubstepping==1){
    380380                        dummy_input = basalelement->GetInput2(SmbRunoffEnum); _assert_(dummy_input);
    381381                }
    382                 _assert_(dummy_input->InstanceEnum()==TransientInput2Enum);
    383                 surface_runoff_input=xDynamicCast<TransientInput2*>(dummy_input); _assert_(surface_runoff_input);
     382                else if(smbsubstepping>1 && smbsubstepping<=hydrologysubstepping){
     383                        dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum, time); _assert_(dummy_input);
     384                }
     385                else{
     386                        dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time); _assert_(dummy_input);
     387                }
     388                surface_runoff_input=xDynamicCast<Input2*>(dummy_input); _assert_(surface_runoff_input);
    384389        }
    385390

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

    r24504 r24565  
    157157                ElementMatrix*     NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
    158158                ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
    159                 void               PicoUpdateBoxid(int* pmax_boxid_basin); 
     159                void               PicoUpdateBoxid(int* pmax_boxid_basin);
    160160                void               PicoUpdateBox(int loopboxid);
    161                 void               PicoComputeBasalMelt(); 
     161                void               PicoComputeBasalMelt();
    162162                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
    163163                void               PositiveDegreeDaySicopolis(bool isfirnwarming);
    … …  
    251251                virtual Input2*    GetInput2(int inputenum)=0;
    252252                virtual Input2*    GetInput2(int inputenum,IssmDouble time)=0;
     253                virtual Input2*    GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time)=0;
    253254                virtual void       GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
    254255                virtual void       GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
    … …  
    340341                virtual Element*   SpawnTopElement(void)=0;
    341342                virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0;
    342                 virtual void       StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa)=0;         
     343                virtual void       StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa)=0;
    343344                virtual void        StrainRateparallel(void)=0;
    344345                virtual void        StrainRateperpendicular(void)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r24481 r24565  
    1 /*! \file Penta.h 
     1/*! \file Penta.h
    22 *  \brief: header file for penta object
    33 */
    … …  
    8585                Input2*        GetInput2(int enumtype);
    8686                Input2*        GetInput2(int enumtype,IssmDouble time);
     87                Input2*        GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
    8788                void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
    8889                void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r24481 r24565  
    1 /*! \file Seg.h 
     1/*! \file Seg.h
    22 *  \brief: header file for seg object
    33 */
    … …  
    6565                Input2*     GetInput2(int enumtype);
    6666                Input2*     GetInput2(int enumtype,IssmDouble time);
     67                Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
    6768                void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
    6869                void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r24481 r24565  
    1 /*! \file Tetra.h 
     1/*! \file Tetra.h
    22 *  \brief: header file for seg object
    33 */
    … …  
    6969                Input2*     GetInput2(int enumtype);
    7070                Input2*     GetInput2(int enumtype,IssmDouble time);
     71                Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
    7172                void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
    7273                void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);

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

    r24505 r24565  
    19111911        }
    19121912}/*}}}*/
     1913Input2*    Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time){/*{{{*/
     1914
     1915        /*Get Input from dataset*/
     1916        if(this->iscollapsed){
     1917                _error_("Get Average input not implemented in Penta yet");
     1918        }
     1919        else{
     1920                TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time);
     1921                if(!input) return input;
     1922
     1923                this->InputServe(input);
     1924                return input;
     1925        }
     1926}/*}}}*/
    19131927void       Tria::GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
    19141928
    … …  
    21022116        IssmDouble dt;
    21032117        int        found,start_offset,end_offset;
    2104         int        averaging_method = 0;
     2118        int        averaging_method=0;
     2119
    21052120
    21062121        /*Get transient input time steps*/
    … …  
    21212136        int offset = start_offset;
    21222137        while(offset < end_offset ){
    2123 
    21242138                if(offset==-1){
    21252139                        /*get values for the first time: */
    … …  
    21322146                                input->GetInputValue(&current_values[iv],gauss);
    21332147                        }
     2148                        dt = timesteps[0] - start_time;
    21342149                }
    21352150                else{
    … …  
    21412156                                input->GetInputValue(&current_values[iv],gauss);
    21422157                        }
    2143                 }
    2144 
    2145                 /*Step between current offset and next*/
    2146                 if(offset==-1){
    2147                         dt = timesteps[0] - start_time;
    2148                 }
    2149                 else if(offset == numtimesteps-1){
    2150                         dt = end_time - timesteps[offset];
    2151                 }
    2152                 else{
    2153                         dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.);
     2158                        if(offset == numtimesteps-1){
     2159                                dt = end_time - timesteps[offset];
     2160                        }
     2161                        else{
     2162                                dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.);
     2163                        }
    21542164                }
    21552165
    … …  
    22202230        IssmDouble *timesteps    = NULL;
    22212231        TransientInput2* transient_input  = this->inputs2->GetTransientInput(input_enum);
     2232
    22222233        transient_input->GetAllTimes(&timesteps,&numtimesteps);
    22232234
    … …  
    55575568        /*Computational flags:*/
    55585569        int bp_compute_fingerprints= 0;
    5559        
     5570
    55605571        /*some paramters first: */
    55615572        this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
    5562        
     5573
    55635574        if(!IsOceanInElement()){
    5564                 /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested 
     5575                /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
    55655576                 *bottom pressure fingerprints:*/
    55665577                if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
    … …  
    56035614        bool computeelastic= true;
    56045615        bool scaleoceanarea= false;
    5605        
     5616
    56065617        /*early return if we are not on an ice cap:*/
    56075618        if(!IsIceOnlyInElement()){
    … …  
    58055816        bool scaleoceanarea= false;
    58065817        bool bp_compute_fingerprints= false;
    5807        
     5818
    58085819        /*we are here to compute fingerprints originating fromn bottom pressure loads:*/
    58095820

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r24505 r24565  
    1 /*! \file Tria.h 
     1/*! \file Tria.h
    22 *  \brief: header file for tria object
    33 */
    … …  
    165165                IssmDouble OceanArea(void);
    166166                IssmDouble OceanAverage(IssmDouble* Sg);
    167                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea); 
     167                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea);
    168168                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    169169                void    SealevelriseEustaticIce(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    … …  
    189189                Input2*        GetInput2(int enumtype);
    190190                Input2*        GetInput2(int enumtype,IssmDouble time);
     191                Input2*        GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time);
    191192                DatasetInput2* GetDatasetInput2(int inputenum);
    192193                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);

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

    r24529 r24565  
    27672767/*}}}*/
    27682768void FemModel::ThicknessAverage(){/*{{{*/
    2769    
     2769
    27702770        int elementswidth                   = this->GetElementsWidth();//just 2D mesh, tria elements
    27712771   int numberofvertices                = this->vertices->NumberOfVertices();//total number of vertices
    … …  
    27802780   for(int i=0;i<this->elements->Size();i++){
    27812781      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
    2782      
     2782
    27832783                /*check if there is ice in this element*/
    2784                 if(!element->IsIceInElement()) continue; 
    2785        
     2784                if(!element->IsIceInElement()) continue;
     2785
    27862786                /*get H on the vertices*/
    27872787                element->GetInputListOnVertices(H,ThicknessEnum);
    … …  
    52275227}
    52285228/*}}}*/
    5229 void FemModel::InitTransientOutputx(int* transientinput_enum,int numoutputs){ /*{{{*/
     5229void FemModel::InitTransientInputx(int* transientinput_enum,int numoutputs){ /*{{{*/
    52305230
    52315231        for(int i=0;i<numoutputs;i++){
    52325232                this->inputs2->DeleteInput(transientinput_enum[i]);
    52335233                this->inputs2->SetTransientInput(transientinput_enum[i],NULL,0);
    5234 
    52355234                /*We need to configure this input!*/
    52365235                TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]); _assert_(transientinput);
    … …  
    52395238}
    52405239/*}}}*/
    5241 void FemModel::StackTransientOutputx(int* input_enum,int* transientinput_enum,IssmDouble subtime,int numoutputs){ /*{{{*/
     5240void FemModel::StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble subtime,int numoutputs){ /*{{{*/
    52425241
    52435242  for(int i=0;i<numoutputs;i++){
    … …  
    52505249                                /*Get the right transient input*/
    52515250                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
    5252                                 TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]); _assert_(transientinput);
     5251                                TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]);
    52535252
    52545253                                /*Get values and lid list*/
    … …  
    52605259
    52615260                                switch(element->ObjectEnum()){
    5262                                         case TriaEnum:  transientinput->AddTriaTimeInput( subtime,numvertices,vertexlids,values,P1Enum); break;
     5261                                        case TriaEnum:  transientinput->AddTriaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); break;
    52635262                                        case PentaEnum: transientinput->AddPentaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); break;
    52645263                                        default: _error_("Not implemented yet");
    … …  
    52715270}
    52725271/*}}}*/
    5273 void FemModel::AverageTransientOutputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs){ /*{{{*/
     5272void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs){ /*{{{*/
    52745273
    52755274        for(int i=0;i<numoutputs;i++){
    … …  
    55955594        int vid,v1,v2,v3;
    55965595        bool refine;
    5597        
     5596
    55985597
    55995598        /*Fill variables*/
    … …  
    56175616        if(!error_elements) _error_("error_elements is NULL!\n");
    56185617        if(groupthreshold<DBL_EPSILON) _error_("group threshold is too small!");
    5619        
     5618
    56205619        /*Get mesh*/
    56215620        this->GetMesh(&index,&x,&y,&numberofvertices,&numberofelements);

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

    r24490 r24565  
    174174                void UpdateConstraintsExtrudeFromTopx();
    175175                void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
    176                 void InitTransientOutputx(int* transientinput_enum,int numoutputs);
    177                 void StackTransientOutputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
    178                 void AverageTransientOutputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs);
     176                void InitTransientInputx(int* transientinput_enum,int numoutputs);
     177                void StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
     178                void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs);
    179179                void UpdateConstraintsx(void);
    180180                int  UpdateVertexPositionsx(void);

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

    r24365 r24565  
    344344        }
    345345}/*}}}*/
     346TriaInput2* Inputs2::GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time){/*{{{*/
     347
     348        /*Get input id*/
     349        int id = EnumToIndex(enum_in);
     350
     351        /*Check that it has the right format*/
     352        Input2* input = this->inputs[id];
     353        if(!input) return NULL;
     354
     355        if(input->ObjectEnum()==TransientInput2Enum){
     356                return xDynamicCast<TransientInput2*>(input)->GetTriaInput(start_time,end_time);
     357        }
     358        else{
     359                _error_("Input "<<EnumToStringx(enum_in)<<" is not an TransientInput2");
     360        }
     361}/*}}}*/
    346362PentaInput2* Inputs2::GetPentaInput(int enum_in){/*{{{*/
    347363
    … …  
    683699
    684700        /*This one only supports P0 and P1 because it assumes col=0*/
    685         _assert_(interpolation==P0Enum || interpolation==P1Enum); 
     701        _assert_(interpolation==P0Enum || interpolation==P1Enum);
    686702
    687703        /*Get input id*/
    … …  
    802818
    803819        /*This one only supports P0 and P1 because it assumes col=0*/
    804         _assert_(interpolation==P0Enum || interpolation==P1Enum); 
     820        _assert_(interpolation==P0Enum || interpolation==P1Enum);
    805821
    806822        /*Get input id*/

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

    r24360 r24565  
    1717#define NUMINPUTS InputsENDEnum - InputsSTARTEnum -1
    1818
    19 /*!\brief Declaration of Inputs class. 
     19/*!\brief Declaration of Inputs class.
    2020 *
    2121 * Declaration of Inputs class.  Inputs are a static array of Input objects.
    22  */ 
     22 */
    2323class Inputs2{
    2424
    … …  
    5454                TriaInput2*      GetTriaInput(int enum_type);
    5555                TriaInput2*      GetTriaInput(int enum_type,IssmDouble time);
     56                TriaInput2*      GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time);
    5657                PentaInput2*     GetPentaInput(int enum_type);
    5758                PentaInput2*     GetPentaInput(int enum_type,IssmDouble time);

  • issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp

    r24351 r24565  
    168168                xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1);
    169169                xDelete(old_timesteps);
     170                xDelete<Input2*>(old_inputs);
    170171        }
    171172
    … …  
    243244        /*Set current time input*/
    244245        this->SetCurrentTimeInput(time);
     246        _assert_(this->current_input);
     247
     248        /*Cast and return*/
     249        if(this->current_input->ObjectEnum()!=TriaInput2Enum){
     250                _error_("Cannot return a TriaInput2");
     251        }
     252        return xDynamicCast<TriaInput2*>(this->current_input);
     253
     254}
     255/*}}}*/
     256TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time){/*{{{*/
     257
     258        /*Set current time input*/
     259        this->SetAverageAsCurrentTimeInput(start_time,end_time);
    245260        _assert_(this->current_input);
    246261
    … …  
    370385
    371386}/*}}}*/
     387void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time){/*{{{*/
     388
     389        IssmDouble dt;
     390        IssmDouble eps=1.0e-6;
     391        IssmDouble dtsum=0;
     392        int        found,start_offset,end_offset;
     393        int        averaging_method = 0;
     394
     395        /*go through the timesteps, and grab offset for start and end*/
     396        found=binary_search(&start_offset,start_time-eps,this->timesteps,this->numtimesteps);
     397        if(!found) _error_("Input not found (is TransientInput sorted ?)");
     398        found=binary_search(&end_offset,end_time+eps,this->timesteps,this->numtimesteps);
     399        if(!found) _error_("Input not found (is TransientInput sorted ?)");
     400
     401        int offset=start_offset;
     402        if(this->current_input) delete this->current_input;
     403        while(offset < end_offset){
     404                if (offset==-1){
     405                        dt=this->timesteps[0]-start_time;
     406                        _assert_(start_time<this->timesteps[0]);
     407                }
     408                else if(offset==this->numtimesteps-1)dt=end_time-this->timesteps[offset];
     409                else if(offset==start_offset && this->timesteps[offset]<start_time) dt=this->timesteps[offset+1]-start_time;
     410                else if(offset==end_offset && this->timesteps[offset]>end_time) dt=end_time-this->timesteps[offset];
     411                else dt=this->timesteps[offset+1]-this->timesteps[offset];
     412                _assert_(dt>0.);
     413
     414                Input2* stepinput=this->inputs[offset+1];
     415
     416                switch(averaging_method){
     417                        case 0: /*Arithmetic mean*/
     418                                if(offset==start_offset){
     419                                        this->current_input = stepinput->copy();
     420                                        this->current_input->Scale(dt);
     421                                }
     422                                else{
     423                                        this->current_input->AXPY(stepinput,dt);
     424                                }
     425                                break;
     426                        case 1: /*Geometric mean*/
     427                                _error_("Geometric not implemented yet");
     428                        case 2: /*Harmonic mean*/
     429                                _error_("Harmonic not implemented yet");
     430                        default:
     431                                _error_("averaging method is not recognised");
     432                }
     433                dtsum+=dt;
     434                offset+=1;
     435        }
     436                /*Integration done, now normalize*/
     437        switch(averaging_method){
     438                case 0: //Arithmetic mean
     439                        this->current_input->Scale(1/(dtsum));
     440                        break;
     441                case 1: /*Geometric mean*/
     442                        _error_("Geometric not implemented yet");
     443                case 2: /*Harmonic mean*/
     444                        _error_("Harmonic not implemented yet");
     445                default:
     446                        _error_("averaging method is not recognised");
     447        }
     448}/*}}}*/
    372449IssmDouble  TransientInput2::GetTimeByOffset(int offset){/*{{{*/
    373450        if(offset<0) offset=0;

  • issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h

    r24351 r24565  
    5050                TriaInput2*  GetTriaInput();
    5151                TriaInput2*  GetTriaInput(IssmDouble time);
     52                TriaInput2*  GetTriaInput(IssmDouble start_time,IssmDouble end_time);
    5253                TriaInput2*  GetTriaInput(int offset);
    5354                PentaInput2* GetPentaInput();
    … …  
    5859                int          GetTimeInputOffset(IssmDouble time);
    5960                void         SetCurrentTimeInput(IssmDouble time);
     61                void         SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time);
    6062                /*numerics:*/
    6163

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

    r24335 r24565  
    9797                                        averagedinput.assign(averagelist,averagelist+2);
    9898                                }
    99                                 femmodel->InitTransientOutputx(&transientinput[0],numaveragedinput);
     99                                femmodel->InitTransientInputx(&transientinput[0],numaveragedinput);
    100100                                while(substep<dtslices){ //loop on hydro dts
    101101                                        substep+=1;
    … …  
    114114                                        solutionsequence_hydro_nonlinear(femmodel);
    115115               /*If we have a sub-timestep we store the substep inputs in a transient input here*/
    116                                         femmodel->StackTransientOutputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
     116                                        femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
    117117                                }
    118118                                /*averaging the stack*/
    119                                 femmodel->AverageTransientOutputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
     119                                femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
    120120
    121121                                /*And reseting to global time*/

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

    r24507 r24565  
    6868                femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum);
    6969
    70                 femmodel->InitTransientOutputx(&transientinput[0],numaveragedinput);
     70                femmodel->InitTransientInputx(&transientinput[0],numaveragedinput);
    7171                while(substep<dtslices){ //loop on sub dts
    7272                        substep+=1;
    … …  
    7979                        analysis->Core(femmodel);
    8080         /*If we have a sub-timestep we store the substep inputs in a transient input here*/
    81          femmodel->StackTransientOutputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
     81         femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
    8282                        delete analysis;
    8383                }
    8484      /*averaging the transient input*/
    85                 femmodel->AverageTransientOutputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
     85                femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
    8686                /*and reset timesteping variables to original*/
    8787                femmodel->parameters->SetParam(global_time,TimeEnum);

  • issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py

    r24385 r24565  
    268268                tensors = [field for field in fieldnames if field[-2:] in ['xx', 'yy', 'xy', 'zz', 'xz', 'yz']]
    269269                non_tensor = [field for field in fieldnames if field not in tensors]
    270                 vectors = [field for field in non_tensor if field[-1] in ['x', 'y', 'z']]
     270                vectors = [field for field in non_tensor if field[-1] in ['x', 'y', 'z'] and field[-4:] not in ['Flux']]
    271271
    272272    #check which field is a real result and print

  • issm/trunk-jpl/src/m/plot/applyoptions.py

    r24256 r24565  
    236236            colorbarfontsize = options.getfieldvalue('colorbarfontsize')
    237237            cb.ax.tick_params(labelsize=colorbarfontsize)
    238     # cb.set_ticks([0, -10])
    239     # cb.set_ticklabels([-10, 0, 10])
    240238        if options.exist('colorbarticks'):
    241239            colorbarticks = options.getfieldvalue('colorbarticks')

  • issm/trunk-jpl/src/m/plot/plot_contour.py

    r24213 r24565  
    3434    norm = options.getfieldvalue('colornorm')
    3535    colors = options.getfieldvalue('contourcolors', 'y')
    36     linestyles = options.getfieldvalue('contourlinestyles', ' - ')
     36    linestyles = options.getfieldvalue('contourlinestyles', '-')
    3737    linewidths = options.getfieldvalue('contourlinewidths', 1)
    3838

  • issm/trunk-jpl/src/m/plot/plot_unit.py

    r24269 r24565  
    6868    # }}}
    6969    # {{{ Get the colormap limits
     70    dataspread = np.nanmax(data) - np.nanmin(data)
     71    if dataspread != 0.:
     72        limextent = np.abs(dataspread / np.nanmean(data))
     73    else:
     74        limextent = 0.
    7075    if options.exist('clim'):
    71         lims = options.getfieldvalue('clim', [np.amin(data), np.amax(data)])
     76        lims = options.getfieldvalue('clim', [np.nanmin(data), np.nanmax(data)])
    7277    elif options.exist('caxis'):
    73         lims = options.getfieldvalue('caxis', [np.amin(data), np.amax(data)])
    74     else:
    75         if np.amin(data) == np.amax(data):
    76             lims = [np.amin(data) * 0.9, np.amax(data) * 1.1]
    77         else:
    78             lims = [np.amin(data), np.amax(data)]
     78        lims = options.getfieldvalue('caxis', [np.nanmin(data), np.nanmax(data)])
     79    else:
     80        if np.nanmin(data) == np.nanmax(data):
     81            lims = [np.nanmin(data) * 0.9, np.nanmax(data) * 1.1]
     82        elif limextent < 1.0e-12:
     83            lims = [np.nanmin(data) - 2 * dataspread, np.nanmax(data) + 2 * dataspread]
     84        else:
     85            lims = [np.nanmin(data), np.nanmax(data)]
    7986    # }}}
    8087    # {{{ Set the spread of the colormap (default is normal
    … …  
    103110            #first deal with colormap
    104111            loccmap = plt.cm.ScalarMappable(cmap=cmap)
    105             loccmap.set_array([np.nanmin(data), np.nanmax(data)])
    106             loccmap.set_clim(vmin=np.nanmin(data), vmax=np.nanmax(data))
     112            loccmap.set_array(lims)
     113            loccmap.set_clim(vmin=lims[0], vmax=lims[1])
    107114
    108115    #dealing with prism sides
    … …  
    173180            #first deal with the colormap
    174181            loccmap = plt.cm.ScalarMappable(cmap=cmap)
    175             loccmap.set_array([np.nanmin(data), np.nanmax(data)])
    176             loccmap.set_clim(vmin=np.nanmin(data), vmax=np.nanmax(data))
     182            loccmap.set_array(lims)
     183            loccmap.set_clim(vmin=lims[0], vmax=lims[1])
    177184
    178185    #deal with prism sides

  • issm/trunk-jpl/test/NightlyRun/runme.py

    r24256 r24565  
    144144                if os.path.isfile(archive_file):
    145145                    os.remove(archive_file)
    146                     for k, fieldname in enumerate(field_names):
    147                         field = np.array(field_values[k], dtype=float)
    148                         if len(field.shape) == 1:
    149                             if np.size(field):
    150                                 field = field.reshape(np.size(field), 1)
    151                             else:
    152                                 field = field.reshape(0, 0)
    153                         elif len(field.shape) == 0:
    154                             field = field.reshape(1, 1)
    155                             # Matlab uses base 1, so use base 1 in labels
    156                         archwrite(archive_file, archive_name + '_field' + str(k + 1), field)
    157                     print(("File {} saved. \n".format(os.path.join('..', 'Archives', archive_name + '.arch'))))
     146                for k, fieldname in enumerate(field_names):
     147                    field = np.array(field_values[k], dtype=float)
     148                    if len(field.shape) == 1:
     149                        if np.size(field):
     150                            field = field.reshape(np.size(field), 1)
     151                        else:
     152                            field = field.reshape(0, 0)
     153                    elif len(field.shape) == 0:
     154                        field = field.reshape(1, 1)
     155                        # Matlab uses base 1, so use base 1 in labels
     156                    archwrite(archive_file, archive_name + '_field' + str(k + 1), field)
     157                print(("File {} saved. \n".format(os.path.join('..', 'Archives', archive_name + '.arch'))))
    158158
    159159                    #ELSE: CHECK TEST
Note: See TracChangeset for help on using the changeset viewer.