Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 24564) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 24565) @@ -335,4 +335,5 @@ int smb_model; int smbsubstepping; + int hydrologysubstepping; IssmDouble dt,scalar,sediment_storing; IssmDouble water_head,sediment_transmitivity; @@ -340,9 +341,9 @@ IssmDouble Jdet; - IssmDouble *xyz_list = NULL; - Input2 *active_element_input = NULL; - Input2 *old_wh_input = NULL; - Input2 *dummy_input = NULL; - TransientInput2 *surface_runoff_input = NULL; + IssmDouble *xyz_list = NULL; + Input2 *active_element_input = NULL; + Input2 *old_wh_input = NULL; + Input2 *dummy_input = NULL; + Input2 *surface_runoff_input = NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -370,16 +371,20 @@ if(dt!= 0.){ - old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input); + old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input); } if(smb_model==SMBgradientscomponentsEnum){ basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum); - if(smbsubstepping>1) { - dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum); _assert_(dummy_input); - } - else { + basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum); + + if(smbsubstepping==1){ dummy_input = basalelement->GetInput2(SmbRunoffEnum); _assert_(dummy_input); } - _assert_(dummy_input->InstanceEnum()==TransientInput2Enum); - surface_runoff_input=xDynamicCast(dummy_input); _assert_(surface_runoff_input); + else if(smbsubstepping>1 && smbsubstepping<=hydrologysubstepping){ + dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum, time); _assert_(dummy_input); + } + else{ + dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time); _assert_(dummy_input); + } + surface_runoff_input=xDynamicCast(dummy_input); _assert_(surface_runoff_input); } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24564) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24565) @@ -157,7 +157,7 @@ ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum); - void PicoUpdateBoxid(int* pmax_boxid_basin); + void PicoUpdateBoxid(int* pmax_boxid_basin); void PicoUpdateBox(int loopboxid); - void PicoComputeBasalMelt(); + void PicoComputeBasalMelt(); void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac); void PositiveDegreeDaySicopolis(bool isfirnwarming); @@ -251,4 +251,5 @@ virtual Input2* GetInput2(int inputenum)=0; virtual Input2* GetInput2(int inputenum,IssmDouble time)=0; + virtual Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time)=0; virtual void GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");}; virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");}; @@ -340,5 +341,5 @@ virtual Element* SpawnTopElement(void)=0; virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0; - virtual void StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa)=0; + virtual void StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa)=0; virtual void StrainRateparallel(void)=0; virtual void StrainRateperpendicular(void)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24564) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24565) @@ -1,3 +1,3 @@ -/*! \file Penta.h +/*! \file Penta.h * \brief: header file for penta object */ @@ -85,4 +85,5 @@ Input2* GetInput2(int enumtype); Input2* GetInput2(int enumtype,IssmDouble time); + Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");}; void GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value); void GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24564) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24565) @@ -1,3 +1,3 @@ -/*! \file Seg.h +/*! \file Seg.h * \brief: header file for seg object */ @@ -65,4 +65,5 @@ Input2* GetInput2(int enumtype); Input2* GetInput2(int enumtype,IssmDouble time); + Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");}; void GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value); void GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value); Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24564) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24565) @@ -1,3 +1,3 @@ -/*! \file Tetra.h +/*! \file Tetra.h * \brief: header file for seg object */ @@ -69,4 +69,5 @@ Input2* GetInput2(int enumtype); Input2* GetInput2(int enumtype,IssmDouble time); + Input2* GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");}; void GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value); void GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24564) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24565) @@ -1911,4 +1911,18 @@ } }/*}}}*/ +Input2* Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time){/*{{{*/ + + /*Get Input from dataset*/ + if(this->iscollapsed){ + _error_("Get Average input not implemented in Penta yet"); + } + else{ + TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time); + if(!input) return input; + + this->InputServe(input); + return input; + } +}/*}}}*/ void Tria::GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/ @@ -2102,5 +2116,6 @@ IssmDouble dt; int found,start_offset,end_offset; - int averaging_method = 0; + int averaging_method=0; + /*Get transient input time steps*/ @@ -2121,5 +2136,4 @@ int offset = start_offset; while(offset < end_offset ){ - if(offset==-1){ /*get values for the first time: */ @@ -2132,4 +2146,5 @@ input->GetInputValue(¤t_values[iv],gauss); } + dt = timesteps[0] - start_time; } else{ @@ -2141,15 +2156,10 @@ input->GetInputValue(¤t_values[iv],gauss); } - } - - /*Step between current offset and next*/ - if(offset==-1){ - dt = timesteps[0] - start_time; - } - else if(offset == numtimesteps-1){ - dt = end_time - timesteps[offset]; - } - else{ - dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.); + if(offset == numtimesteps-1){ + dt = end_time - timesteps[offset]; + } + else{ + dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.); + } } @@ -2220,4 +2230,5 @@ IssmDouble *timesteps = NULL; TransientInput2* transient_input = this->inputs2->GetTransientInput(input_enum); + transient_input->GetAllTimes(×teps,&numtimesteps); @@ -5557,10 +5568,10 @@ /*Computational flags:*/ int bp_compute_fingerprints= 0; - + /*some paramters first: */ this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); - + if(!IsOceanInElement()){ - /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested + /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested *bottom pressure fingerprints:*/ if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); @@ -5603,5 +5614,5 @@ bool computeelastic= true; bool scaleoceanarea= false; - + /*early return if we are not on an ice cap:*/ if(!IsIceOnlyInElement()){ @@ -5805,5 +5816,5 @@ bool scaleoceanarea= false; bool bp_compute_fingerprints= false; - + /*we are here to compute fingerprints originating fromn bottom pressure loads:*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24564) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24565) @@ -1,3 +1,3 @@ -/*! \file Tria.h +/*! \file Tria.h * \brief: header file for tria object */ @@ -165,5 +165,5 @@ IssmDouble OceanArea(void); IssmDouble OceanAverage(IssmDouble* Sg); - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea); + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea); void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); void SealevelriseEustaticIce(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); @@ -189,4 +189,5 @@ Input2* GetInput2(int enumtype); Input2* GetInput2(int enumtype,IssmDouble time); + Input2* GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time); DatasetInput2* GetDatasetInput2(int inputenum); void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24564) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24565) @@ -2767,5 +2767,5 @@ /*}}}*/ void FemModel::ThicknessAverage(){/*{{{*/ - + int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements int numberofvertices = this->vertices->NumberOfVertices();//total number of vertices @@ -2780,8 +2780,8 @@ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - + /*check if there is ice in this element*/ - if(!element->IsIceInElement()) continue; - + if(!element->IsIceInElement()) continue; + /*get H on the vertices*/ element->GetInputListOnVertices(H,ThicknessEnum); @@ -5227,10 +5227,9 @@ } /*}}}*/ -void FemModel::InitTransientOutputx(int* transientinput_enum,int numoutputs){ /*{{{*/ +void FemModel::InitTransientInputx(int* transientinput_enum,int numoutputs){ /*{{{*/ for(int i=0;iinputs2->DeleteInput(transientinput_enum[i]); this->inputs2->SetTransientInput(transientinput_enum[i],NULL,0); - /*We need to configure this input!*/ TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]); _assert_(transientinput); @@ -5239,5 +5238,5 @@ } /*}}}*/ -void FemModel::StackTransientOutputx(int* input_enum,int* transientinput_enum,IssmDouble subtime,int numoutputs){ /*{{{*/ +void FemModel::StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble subtime,int numoutputs){ /*{{{*/ for(int i=0;i(elements->GetObjectByOffset(j)); - TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]); _assert_(transientinput); + TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]); /*Get values and lid list*/ @@ -5260,5 +5259,5 @@ switch(element->ObjectEnum()){ - case TriaEnum: transientinput->AddTriaTimeInput( subtime,numvertices,vertexlids,values,P1Enum); break; + case TriaEnum: transientinput->AddTriaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); break; case PentaEnum: transientinput->AddPentaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); break; default: _error_("Not implemented yet"); @@ -5271,5 +5270,5 @@ } /*}}}*/ -void FemModel::AverageTransientOutputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs){ /*{{{*/ +void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs){ /*{{{*/ for(int i=0;iGetMesh(&index,&x,&y,&numberofvertices,&numberofelements); Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24564) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24565) @@ -174,7 +174,7 @@ void UpdateConstraintsExtrudeFromTopx(); void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count); - void InitTransientOutputx(int* transientinput_enum,int numoutputs); - void StackTransientOutputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs); - void AverageTransientOutputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs); + void InitTransientInputx(int* transientinput_enum,int numoutputs); + void StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs); + void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs); void UpdateConstraintsx(void); int UpdateVertexPositionsx(void); Index: /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp (revision 24564) +++ /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp (revision 24565) @@ -344,4 +344,20 @@ } }/*}}}*/ +TriaInput2* Inputs2::GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time){/*{{{*/ + + /*Get input id*/ + int id = EnumToIndex(enum_in); + + /*Check that it has the right format*/ + Input2* input = this->inputs[id]; + if(!input) return NULL; + + if(input->ObjectEnum()==TransientInput2Enum){ + return xDynamicCast(input)->GetTriaInput(start_time,end_time); + } + else{ + _error_("Input "<timesteps,old_timesteps,this->numtimesteps-1); xDelete(old_timesteps); + xDelete(old_inputs); } @@ -243,4 +244,18 @@ /*Set current time input*/ this->SetCurrentTimeInput(time); + _assert_(this->current_input); + + /*Cast and return*/ + if(this->current_input->ObjectEnum()!=TriaInput2Enum){ + _error_("Cannot return a TriaInput2"); + } + return xDynamicCast(this->current_input); + +} +/*}}}*/ +TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time){/*{{{*/ + + /*Set current time input*/ + this->SetAverageAsCurrentTimeInput(start_time,end_time); _assert_(this->current_input); @@ -370,4 +385,66 @@ }/*}}}*/ +void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time){/*{{{*/ + + IssmDouble dt; + IssmDouble eps=1.0e-6; + IssmDouble dtsum=0; + int found,start_offset,end_offset; + int averaging_method = 0; + + /*go through the timesteps, and grab offset for start and end*/ + found=binary_search(&start_offset,start_time-eps,this->timesteps,this->numtimesteps); + if(!found) _error_("Input not found (is TransientInput sorted ?)"); + found=binary_search(&end_offset,end_time+eps,this->timesteps,this->numtimesteps); + if(!found) _error_("Input not found (is TransientInput sorted ?)"); + + int offset=start_offset; + if(this->current_input) delete this->current_input; + while(offset < end_offset){ + if (offset==-1){ + dt=this->timesteps[0]-start_time; + _assert_(start_timetimesteps[0]); + } + else if(offset==this->numtimesteps-1)dt=end_time-this->timesteps[offset]; + else if(offset==start_offset && this->timesteps[offset]timesteps[offset+1]-start_time; + else if(offset==end_offset && this->timesteps[offset]>end_time) dt=end_time-this->timesteps[offset]; + else dt=this->timesteps[offset+1]-this->timesteps[offset]; + _assert_(dt>0.); + + Input2* stepinput=this->inputs[offset+1]; + + switch(averaging_method){ + case 0: /*Arithmetic mean*/ + if(offset==start_offset){ + this->current_input = stepinput->copy(); + this->current_input->Scale(dt); + } + else{ + this->current_input->AXPY(stepinput,dt); + } + break; + case 1: /*Geometric mean*/ + _error_("Geometric not implemented yet"); + case 2: /*Harmonic mean*/ + _error_("Harmonic not implemented yet"); + default: + _error_("averaging method is not recognised"); + } + dtsum+=dt; + offset+=1; + } + /*Integration done, now normalize*/ + switch(averaging_method){ + case 0: //Arithmetic mean + this->current_input->Scale(1/(dtsum)); + break; + case 1: /*Geometric mean*/ + _error_("Geometric not implemented yet"); + case 2: /*Harmonic mean*/ + _error_("Harmonic not implemented yet"); + default: + _error_("averaging method is not recognised"); + } +}/*}}}*/ IssmDouble TransientInput2::GetTimeByOffset(int offset){/*{{{*/ if(offset<0) offset=0; Index: /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h (revision 24564) +++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h (revision 24565) @@ -50,4 +50,5 @@ TriaInput2* GetTriaInput(); TriaInput2* GetTriaInput(IssmDouble time); + TriaInput2* GetTriaInput(IssmDouble start_time,IssmDouble end_time); TriaInput2* GetTriaInput(int offset); PentaInput2* GetPentaInput(); @@ -58,4 +59,5 @@ int GetTimeInputOffset(IssmDouble time); void SetCurrentTimeInput(IssmDouble time); + void SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time); /*numerics:*/ Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 24564) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 24565) @@ -97,5 +97,5 @@ averagedinput.assign(averagelist,averagelist+2); } - femmodel->InitTransientOutputx(&transientinput[0],numaveragedinput); + femmodel->InitTransientInputx(&transientinput[0],numaveragedinput); while(substepStackTransientOutputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput); + femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput); } /*averaging the stack*/ - femmodel->AverageTransientOutputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput); + femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput); /*And reseting to global time*/ Index: /issm/trunk-jpl/src/c/cores/smb_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/smb_core.cpp (revision 24564) +++ /issm/trunk-jpl/src/c/cores/smb_core.cpp (revision 24565) @@ -68,5 +68,5 @@ femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum); - femmodel->InitTransientOutputx(&transientinput[0],numaveragedinput); + femmodel->InitTransientInputx(&transientinput[0],numaveragedinput); while(substepCore(femmodel); /*If we have a sub-timestep we store the substep inputs in a transient input here*/ - femmodel->StackTransientOutputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput); + femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput); delete analysis; } /*averaging the transient input*/ - femmodel->AverageTransientOutputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput); + femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput); /*and reset timesteping variables to original*/ femmodel->parameters->SetParam(global_time,TimeEnum); Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 24564) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 24565) @@ -268,5 +268,5 @@ tensors = [field for field in fieldnames if field[-2:] in ['xx', 'yy', 'xy', 'zz', 'xz', 'yz']] non_tensor = [field for field in fieldnames if field not in tensors] - vectors = [field for field in non_tensor if field[-1] in ['x', 'y', 'z']] + vectors = [field for field in non_tensor if field[-1] in ['x', 'y', 'z'] and field[-4:] not in ['Flux']] #check which field is a real result and print Index: /issm/trunk-jpl/src/m/plot/applyoptions.py =================================================================== --- /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 24564) +++ /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 24565) @@ -236,6 +236,4 @@ colorbarfontsize = options.getfieldvalue('colorbarfontsize') cb.ax.tick_params(labelsize=colorbarfontsize) - # cb.set_ticks([0, -10]) - # cb.set_ticklabels([-10, 0, 10]) if options.exist('colorbarticks'): colorbarticks = options.getfieldvalue('colorbarticks') Index: /issm/trunk-jpl/src/m/plot/plot_contour.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_contour.py (revision 24564) +++ /issm/trunk-jpl/src/m/plot/plot_contour.py (revision 24565) @@ -34,5 +34,5 @@ norm = options.getfieldvalue('colornorm') colors = options.getfieldvalue('contourcolors', 'y') - linestyles = options.getfieldvalue('contourlinestyles', ' - ') + linestyles = options.getfieldvalue('contourlinestyles', '-') linewidths = options.getfieldvalue('contourlinewidths', 1) Index: /issm/trunk-jpl/src/m/plot/plot_unit.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 24564) +++ /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 24565) @@ -68,13 +68,20 @@ # }}} # {{{ Get the colormap limits + dataspread = np.nanmax(data) - np.nanmin(data) + if dataspread != 0.: + limextent = np.abs(dataspread / np.nanmean(data)) + else: + limextent = 0. if options.exist('clim'): - lims = options.getfieldvalue('clim', [np.amin(data), np.amax(data)]) + lims = options.getfieldvalue('clim', [np.nanmin(data), np.nanmax(data)]) elif options.exist('caxis'): - lims = options.getfieldvalue('caxis', [np.amin(data), np.amax(data)]) - else: - if np.amin(data) == np.amax(data): - lims = [np.amin(data) * 0.9, np.amax(data) * 1.1] - else: - lims = [np.amin(data), np.amax(data)] + lims = options.getfieldvalue('caxis', [np.nanmin(data), np.nanmax(data)]) + else: + if np.nanmin(data) == np.nanmax(data): + lims = [np.nanmin(data) * 0.9, np.nanmax(data) * 1.1] + elif limextent < 1.0e-12: + lims = [np.nanmin(data) - 2 * dataspread, np.nanmax(data) + 2 * dataspread] + else: + lims = [np.nanmin(data), np.nanmax(data)] # }}} # {{{ Set the spread of the colormap (default is normal @@ -103,6 +110,6 @@ #first deal with colormap loccmap = plt.cm.ScalarMappable(cmap=cmap) - loccmap.set_array([np.nanmin(data), np.nanmax(data)]) - loccmap.set_clim(vmin=np.nanmin(data), vmax=np.nanmax(data)) + loccmap.set_array(lims) + loccmap.set_clim(vmin=lims[0], vmax=lims[1]) #dealing with prism sides @@ -173,6 +180,6 @@ #first deal with the colormap loccmap = plt.cm.ScalarMappable(cmap=cmap) - loccmap.set_array([np.nanmin(data), np.nanmax(data)]) - loccmap.set_clim(vmin=np.nanmin(data), vmax=np.nanmax(data)) + loccmap.set_array(lims) + loccmap.set_clim(vmin=lims[0], vmax=lims[1]) #deal with prism sides Index: /issm/trunk-jpl/test/NightlyRun/runme.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/runme.py (revision 24564) +++ /issm/trunk-jpl/test/NightlyRun/runme.py (revision 24565) @@ -144,16 +144,16 @@ if os.path.isfile(archive_file): os.remove(archive_file) - for k, fieldname in enumerate(field_names): - field = np.array(field_values[k], dtype=float) - if len(field.shape) == 1: - if np.size(field): - field = field.reshape(np.size(field), 1) - else: - field = field.reshape(0, 0) - elif len(field.shape) == 0: - field = field.reshape(1, 1) - # Matlab uses base 1, so use base 1 in labels - archwrite(archive_file, archive_name + '_field' + str(k + 1), field) - print(("File {} saved. \n".format(os.path.join('..', 'Archives', archive_name + '.arch')))) + for k, fieldname in enumerate(field_names): + field = np.array(field_values[k], dtype=float) + if len(field.shape) == 1: + if np.size(field): + field = field.reshape(np.size(field), 1) + else: + field = field.reshape(0, 0) + elif len(field.shape) == 0: + field = field.reshape(1, 1) + # Matlab uses base 1, so use base 1 in labels + archwrite(archive_file, archive_name + '_field' + str(k + 1), field) + print(("File {} saved. \n".format(os.path.join('..', 'Archives', archive_name + '.arch')))) #ELSE: CHECK TEST