Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25757) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25758) @@ -175,7 +175,5 @@ IssmDouble planetradius=0; IssmDouble planetarea=0; - bool rigid=false; bool elastic=false; - bool rotation=false; int numoutputs; @@ -243,5 +241,4 @@ } - /*Deal with dsl multi-model ensembles: {{{*/ iomodel->FetchData(&dslmodel,"md.dsl.model"); @@ -263,16 +260,7 @@ } /*}}}*/ /*Deal with elasticity {{{*/ - iomodel->FetchData(&rigid,"md.solidearth.settings.rigid"); iomodel->FetchData(&elastic,"md.solidearth.settings.elastic"); - iomodel->FetchData(&rotation,"md.solidearth.settings.rotation"); - - if(elastic | rigid){ - /*compute green functions for a range of angles*/ - iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); - M=reCast(180./degacc+1.); - } - - /*love numbers: */ if(elastic){ + /*love numbers: */ iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); @@ -281,4 +269,5 @@ iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); @@ -289,63 +278,44 @@ parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); + /*compute elastic green function for a range of angles*/ + iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); + M=reCast(180./degacc+1.); + // AD performance is sensitive to calls to ensurecontiguous. // // Providing "t" will cause ensurecontiguous to be called. #ifdef _HAVE_AD_ + G_rigid=xNew(M,"t"); G_elastic=xNew(M,"t"); U_elastic=xNew(M,"t"); H_elastic=xNew(M,"t"); #else + G_rigid=xNew(M); G_elastic=xNew(M); U_elastic=xNew(M); H_elastic=xNew(M); #endif - } - if(rigid){ - #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); - #else - G_rigid=xNew(M); - #endif - } - - if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); - - if(rigid | elastic){ /*compute combined legendre + love number (elastic green function:*/ m=DetermineLocalSize(M,IssmComm::GetComm()); GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - } - if(elastic){ + // AD performance is sensitive to calls to ensurecontiguous. + // // Providing "t" will cause ensurecontiguous to be called. #ifdef _HAVE_AD_ G_elastic_local=xNew(m,"t"); + G_rigid_local=xNew(m,"t"); U_elastic_local=xNew(m,"t"); H_elastic_local=xNew(m,"t"); #else G_elastic_local=xNew(m); + G_rigid_local=xNew(m); U_elastic_local=xNew(m); H_elastic_local=xNew(m); #endif - } - if(rigid){ - #ifdef _HAVE_AD_ - G_rigid_local=xNew(m,"t"); - #else - G_rigid_local=xNew(m); - #endif - } - - if(rigid){ + for(int i=lower_row;i(i)*degacc * PI / 180.0; + G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); - } - } - if(elastic){ - for(int i=lower_row;i(i)*degacc * PI / 180.0; - G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; @@ -386,6 +356,4 @@ } } - } - if(rigid){ /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/ @@ -401,50 +369,43 @@ /*All gather:*/ ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - if(elastic){ - ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - } - + ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + /*free resources: */ xDelete(recvcounts); xDelete(displs); + /*}}}*/ /*Avoid singularity at 0: */ G_rigid[0]=G_rigid[1]; + parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); G_elastic[0]=G_elastic[1]; + parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); U_elastic[0]=U_elastic[1]; + parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); H_elastic[0]=H_elastic[1]; - - /*Save our precomputed tables into parameters*/ - parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); - if(elastic){ - parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); - parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); - parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); - } + parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); /*free resources: */ + xDelete(love_h); + xDelete(love_k); + xDelete(love_l); + xDelete(love_th); + xDelete(love_tk); + xDelete(love_tl); xDelete(G_rigid); xDelete(G_rigid_local); - if(elastic){ - xDelete(love_h); - xDelete(love_k); - xDelete(love_l); - xDelete(love_th); - xDelete(love_tk); - xDelete(love_tl); - xDelete(G_elastic); - xDelete(G_elastic_local); - xDelete(U_elastic); - xDelete(U_elastic_local); - xDelete(H_elastic); - xDelete(H_elastic_local); - } + xDelete(G_elastic); + xDelete(G_elastic_local); + xDelete(U_elastic); + xDelete(U_elastic_local); + xDelete(H_elastic); + xDelete(H_elastic_local); } - /*}}}*/ /*Indicate we have not yet run the Geometry Core module: */ parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false)); + /*}}}*/ /*Transitions:{{{ */ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25757) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25758) @@ -5446,6 +5446,6 @@ /*Computational flags:*/ - bool computerigid = false; - bool computeelastic = false; + bool computerigid = true; + bool computeelastic = true; int horiz; @@ -5463,5 +5463,5 @@ this->inputs->SetDoubleInput(AreaEnum,this->lid,area); this->AddInput(SealevelAreaEnum,&area,P0Enum); - if(!computerigid)return; + if(!computerigid & !computeelastic)return; /*recover precomputed green function kernels:*/ @@ -5469,24 +5469,20 @@ parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); + parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + /*allocate indices:*/ indices=xNew(gsize); G=xNewZeroInit(gsize); - - if(computeelastic){ - parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); - - /*allocate indices:*/ - GU=xNewZeroInit(gsize); - if(horiz){ - GN=xNewZeroInit(gsize); - GE=xNewZeroInit(gsize); - } + GU=xNewZeroInit(gsize); + if(horiz){ + GN=xNewZeroInit(gsize); + GE=xNewZeroInit(gsize); } @@ -5503,8 +5499,8 @@ /*compute gravity center in lat,long: */ late= asin(z_element/planetradius); - longe = atan2(y_element,x_element); - /*}}}*/ + longe = atan2(y_element,x_element); constant=3/rho_earth*area/planetarea; + for(int i=0;iinputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); - this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); - if(computeelastic){ - this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); - if(horiz){ - this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); - this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); - } - } - this->inputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); + this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); + this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); + this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); + if(horiz){ + this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); + this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); + } /*Free allocations:*/ #ifdef _HAVE_RESTRICT_ + delete indices; delete G; - if(computeelastic){ - delete GU; - if(horiz){ - delete GN; - delete GE; - } + delete GU; + if(horiz){ + delete GN; + delete GE; } #else + xDelete(indices); xDelete(G); - if(computeelastic){ - xDelete(GU); - if(horiz){ - xDelete(GN); - xDelete(GE); - } + xDelete(GU); + if(horiz){ + xDelete(GN); + xDelete(GE); } #endif @@ -5591,5 +5583,5 @@ bool notfullygrounded=false; bool scaleoceanarea= false; - bool computerigid= false; + bool computeelastic= true; int glfraction=1; @@ -5637,9 +5629,10 @@ rho_ice=FindParam(MaterialsRhoIceEnum); rho_water=FindParam(MaterialsRhoSeawaterEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); + this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); /*retrieve precomputed G:*/ - if(computerigid)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); + if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ @@ -5695,5 +5688,5 @@ _assert_(!xIsNan(bslrice)); - if(computerigid){ + if(computeelastic){ /*convert from m to kg/m^2:*/ I=I*rho_ice*phi; @@ -5721,5 +5714,5 @@ bool notfullygrounded=false; bool scaleoceanarea= false; - bool computerigid= false; + bool computeelastic= true; /*elastic green function:*/ @@ -5751,9 +5744,9 @@ rho_water=FindParam(MaterialsRhoSeawaterEnum); rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); /*retrieve precomputed G:*/ - if(computerigid)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); + if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ @@ -5771,5 +5764,5 @@ _assert_(!xIsNan(bslrhydro)); - if(computerigid){ + if(computeelastic){ /*convert from m to kg/m^2:*/ W=W*rho_freshwater*phi; @@ -5795,5 +5788,5 @@ IssmDouble BP; //change in bottom pressure (Farrel and Clarke, Equ. 4) IssmDouble constant; - bool computeelastic= false; + bool computeelastic= true; /*elastic green function:*/ @@ -5897,5 +5890,5 @@ /*computational flags:*/ - bool computeelastic= false; + bool computeelastic= true; /*early return if we are not on the ocean or on an ice cap:*/ Index: /issm/trunk-jpl/src/m/classes/clusters/generic.py =================================================================== --- /issm/trunk-jpl/src/m/classes/clusters/generic.py (revision 25757) +++ /issm/trunk-jpl/src/m/classes/clusters/generic.py (revision 25758) @@ -1,25 +1,26 @@ -import numpy as np from socket import gethostname from subprocess import call -from IssmConfig import IssmConfig -from issmdir import issmdir -from pairoptions import pairoptions -from issmssh import issmssh -from issmscpin import issmscpin -from issmscpout import issmscpout -from MatlabFuncs import ispc + +import numpy as np + try: from generic_settings import generic_settings except ImportError: print('Warning generic_settings.py not found, default will be used') +from MatlabFuncs import ispc +from IssmConfig import IssmConfig +from issmdir import issmdir +from issmscpin import issmscpin +from issmscpout import issmscpout +from issmssh import issmssh +from pairoptions import pairoptions class generic(object): - """ - GENERIC cluster class definition - - Usage: - cluster = generic('name', 'astrid', 'np', 3) - cluster = generic('name', gethostname(), 'np', 3, 'login', 'username') + """GENERIC cluster class definition + + Usage: + cluster = generic('name', 'astrid', 'np', 3) + cluster = generic('name', gethostname(), 'np', 3, 'login', 'username') """ @@ -50,6 +51,4 @@ except NameError: print("generic_settings.py not found, using default settings") - # else: - # raise #OK get other fields Index: /issm/trunk-jpl/src/m/classes/mask.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mask.m (revision 25757) +++ /issm/trunk-jpl/src/m/classes/mask.m (revision 25758) @@ -114,5 +114,5 @@ WriteData(fid,prefix,'object',self,'fieldname','ice_levelset','name','md.mask.ice_levelset','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); else - WriteData(fid,prefix,'object',self,'fieldname','ice_levelset','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','ice_levelset','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); end end % }}} Index: /issm/trunk-jpl/src/m/classes/mask.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mask.py (revision 25757) +++ /issm/trunk-jpl/src/m/classes/mask.py (revision 25758) @@ -8,10 +8,10 @@ class mask(object): - ''' - MASK class definition + """MASK class definition - Usage: - mask = mask() - ''' + Usage: + mask = mask() + """ + def __init__(self): # {{{ self.ice_levelset = float('NaN') @@ -91,5 +91,5 @@ def marshall(self, prefix, md, fid): # {{{ - WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_levelset', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid, prefix, 'object', self, 'fieldname', 'ice_levelset', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_levelset', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'ice_levelset', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) # }}} Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 25757) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 25758) @@ -192,14 +192,15 @@ methods function md = model(varargin) % {{{ - switch nargin case 0 - md=setdefaultparameters(md,'earth'); + md=setdefaultparameters(md); + case 1 + error('model constructor not supported yet'); + otherwise - options=pairoptions(varargin{:}); - planet=getfieldvalue(options,'planet','earth'); - md=setdefaultparameters(md,planet); - end + error('model constructor error message: 0 of 1 argument only in input.'); + end + end %}}} @@ -424,5 +425,5 @@ if isobject(md.outputdefinition.definitions{i}) %get subfields - solutionsubfields=fieldnames(md.outputdefinition.definitions{i}); + solutionsubfields=fields(md.outputdefinition.definitions{i}); for j=1:length(solutionsubfields), field=md.outputdefinition.definitions{i}.(solutionsubfields{j}); @@ -548,5 +549,5 @@ %loop over model fields - model_fields=fieldnames(md1); + model_fields=fields(md1); for i=1:length(model_fields), %get field @@ -554,5 +555,5 @@ fieldsize=size(field); if isobject(field), %recursive call - object_fields=fieldnames(md1.(model_fields{i})); + object_fields=fields(md1.(model_fields{i})); for j=1:length(object_fields), %get field @@ -721,5 +722,5 @@ if isstruct(md1.results), md2.results=struct(); - solutionfields=fieldnames(md1.results); + solutionfields=fields(md1.results); for i=1:length(solutionfields), if isstruct(md1.results.(solutionfields{i})) @@ -728,5 +729,5 @@ for p=1:length(md1.results.(solutionfields{i})) current = md1.results.(solutionfields{i})(p); - solutionsubfields=fieldnames(current); + solutionsubfields=fields(current); for j=1:length(solutionsubfields), field=md1.results.(solutionfields{i})(p).(solutionsubfields{j}); @@ -757,5 +758,5 @@ if isobject(md1.outputdefinition.definitions{i}) %get subfields - solutionsubfields=fieldnames(md1.outputdefinition.definitions{i}); + solutionsubfields=fields(md1.outputdefinition.definitions{i}); for j=1:length(solutionsubfields), field=md1.outputdefinition.definitions{i}.(solutionsubfields{j}); @@ -1287,5 +1288,5 @@ end end% }}} - function md = setdefaultparameters(md,planet) % {{{ + function md = setdefaultparameters(md) % {{{ %initialize subclasses @@ -1299,5 +1300,5 @@ md.friction = friction(); md.rifts = rifts(); - md.solidearth = solidearth(planet); + md.solidearth = solidearth(); md.dsl = dsl(); md.timestepping = timestepping(); Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 25757) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 25758) @@ -173,48 +173,53 @@ def __repr__(obj): #{{{ + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # + #print "Here %s the number: %d" % ("is", 37) - string = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties") - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("geometry", "[%s %s]" % ("1x1", obj.geometry.__class__.__name__), "surface elevation, bedrock topography, ice thickness, ...")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("constants", "[%s %s]" % ("1x1", obj.constants.__class__.__name__), "physical constants")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("smb", "[%s %s]" % ("1x1", obj.smb.__class__.__name__), "surface mass balance")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("basalforcings", "[%s %s]" % ("1x1", obj.basalforcings.__class__.__name__), "bed forcings")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("materials", "[%s %s]" % ("1x1", obj.materials.__class__.__name__), "material properties")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("damage", "[%s %s]" % ("1x1", obj.damage.__class__.__name__), "damage propagation laws")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("friction", "[%s %s]" % ("1x1", obj.friction.__class__.__name__), "basal friction / drag properties")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("flowequation", "[%s %s]" % ("1x1", obj.flowequation.__class__.__name__), "flow equations")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("timestepping", "[%s %s]" % ("1x1", obj.timestepping.__class__.__name__), "time stepping for transient models")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("initialization", "[%s %s]" % ("1x1", obj.initialization.__class__.__name__), "initial guess / state")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("rifts", "[%s %s]" % ("1x1", obj.rifts.__class__.__name__), "rifts properties")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("solidearth", "[%s %s]" % ("1x1", obj.solidearth.__class__.__name__), "solidearth inputs and settings")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("dsl", "[%s %s]" % ("1x1", obj.dsl.__class__.__name__), "dynamic sea level")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("debug", "[%s %s]" % ("1x1", obj.debug.__class__.__name__), "debugging tools (valgrind, gprof)")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("verbose", "[%s %s]" % ("1x1", obj.verbose.__class__.__name__), "verbosity level in solve")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("settings", "[%s %s]" % ("1x1", obj.settings.__class__.__name__), "settings properties")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("toolkits", "[%s %s]" % ("1x1", obj.toolkits.__class__.__name__), "PETSc options for each solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("cluster", "[%s %s]" % ("1x1", obj.cluster.__class__.__name__), "cluster parameters (number of cpus...)")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("balancethickness", "[%s %s]" % ("1x1", obj.balancethickness.__class__.__name__), "parameters for balancethickness solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("stressbalance", "[%s %s]" % ("1x1", obj.stressbalance.__class__.__name__), "parameters for stressbalance solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("groundingline", "[%s %s]" % ("1x1", obj.groundingline.__class__.__name__), "parameters for groundingline solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("hydrology", "[%s %s]" % ("1x1", obj.hydrology.__class__.__name__), "parameters for hydrology solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("masstransport", "[%s %s]" % ("1x1", obj.masstransport.__class__.__name__), "parameters for masstransport solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("thermal", "[%s %s]" % ("1x1", obj.thermal.__class__.__name__), "parameters for thermal solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("steadystate", "[%s %s]" % ("1x1", obj.steadystate.__class__.__name__), "parameters for steadystate solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("transient", "[%s %s]" % ("1x1", obj.transient.__class__.__name__), "parameters for transient solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("levelset", "[%s %s]" % ("1x1", obj.levelset.__class__.__name__), "parameters for moving boundaries (level - set method)")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution")) - string = "%s\n%s" % (string, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution")) - string = "%s\n%s" % (string, '%19s: %-22s -- %s' % ("love", "[%s %s]" % ("1x1", obj.love.__class__.__name__), "parameters for love solution")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("autodiff", "[%s %s]" % ("1x1", obj.autodiff.__class__.__name__), "automatic differentiation parameters")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("inversion", "[%s %s]" % ("1x1", obj.inversion.__class__.__name__), "parameters for inverse methods")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("qmu", "[%s %s]" % ("1x1", obj.qmu.__class__.__name__), "dakota properties")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("amr", "[%s %s]" % ("1x1", obj.amr.__class__.__name__), "adaptive mesh refinement properties")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("outputdefinition", "[%s %s]" % ("1x1", obj.outputdefinition.__class__.__name__), "output definition")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("results", "[%s %s]" % ("1x1", obj.results.__class__.__name__), "model results")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("radaroverlay", "[%s %s]" % ("1x1", obj.radaroverlay.__class__.__name__), "radar image for plot overlay")) - string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields")) - return string + s = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties") + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("geometry", "[%s %s]" % ("1x1", obj.geometry.__class__.__name__), "surface elevation, bedrock topography, ice thickness, ...")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("constants", "[%s %s]" % ("1x1", obj.constants.__class__.__name__), "physical constants")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("smb", "[%s %s]" % ("1x1", obj.smb.__class__.__name__), "surface mass balance")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("basalforcings", "[%s %s]" % ("1x1", obj.basalforcings.__class__.__name__), "bed forcings")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("materials", "[%s %s]" % ("1x1", obj.materials.__class__.__name__), "material properties")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("damage", "[%s %s]" % ("1x1", obj.damage.__class__.__name__), "damage propagation laws")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("friction", "[%s %s]" % ("1x1", obj.friction.__class__.__name__), "basal friction / drag properties")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("flowequation", "[%s %s]" % ("1x1", obj.flowequation.__class__.__name__), "flow equations")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("timestepping", "[%s %s]" % ("1x1", obj.timestepping.__class__.__name__), "time stepping for transient models")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("initialization", "[%s %s]" % ("1x1", obj.initialization.__class__.__name__), "initial guess / state")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("rifts", "[%s %s]" % ("1x1", obj.rifts.__class__.__name__), "rifts properties")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("solidearth", "[%s %s]" % ("1x1", obj.solidearth.__class__.__name__), "solidearth inputs and settings")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("dsl", "[%s %s]" % ("1x1", obj.dsl.__class__.__name__), "dynamic sea level")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("debug", "[%s %s]" % ("1x1", obj.debug.__class__.__name__), "debugging tools (valgrind, gprof)")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("verbose", "[%s %s]" % ("1x1", obj.verbose.__class__.__name__), "verbosity level in solve")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("settings", "[%s %s]" % ("1x1", obj.settings.__class__.__name__), "settings properties")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("toolkits", "[%s %s]" % ("1x1", obj.toolkits.__class__.__name__), "PETSc options for each solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("cluster", "[%s %s]" % ("1x1", obj.cluster.__class__.__name__), "cluster parameters (number of cpus...)")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("balancethickness", "[%s %s]" % ("1x1", obj.balancethickness.__class__.__name__), "parameters for balancethickness solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("stressbalance", "[%s %s]" % ("1x1", obj.stressbalance.__class__.__name__), "parameters for stressbalance solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("groundingline", "[%s %s]" % ("1x1", obj.groundingline.__class__.__name__), "parameters for groundingline solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("hydrology", "[%s %s]" % ("1x1", obj.hydrology.__class__.__name__), "parameters for hydrology solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("masstransport", "[%s %s]" % ("1x1", obj.masstransport.__class__.__name__), "parameters for masstransport solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("thermal", "[%s %s]" % ("1x1", obj.thermal.__class__.__name__), "parameters for thermal solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("steadystate", "[%s %s]" % ("1x1", obj.steadystate.__class__.__name__), "parameters for steadystate solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("transient", "[%s %s]" % ("1x1", obj.transient.__class__.__name__), "parameters for transient solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("levelset", "[%s %s]" % ("1x1", obj.levelset.__class__.__name__), "parameters for moving boundaries (level - set method)")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution")) + s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution")) + s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("love", "[%s %s]" % ("1x1", obj.love.__class__.__name__), "parameters for love solution")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("autodiff", "[%s %s]" % ("1x1", obj.autodiff.__class__.__name__), "automatic differentiation parameters")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("inversion", "[%s %s]" % ("1x1", obj.inversion.__class__.__name__), "parameters for inverse methods")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("qmu", "[%s %s]" % ("1x1", obj.qmu.__class__.__name__), "dakota properties")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("amr", "[%s %s]" % ("1x1", obj.amr.__class__.__name__), "adaptive mesh refinement properties")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("outputdefinition", "[%s %s]" % ("1x1", obj.outputdefinition.__class__.__name__), "output definition")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("results", "[%s %s]" % ("1x1", obj.results.__class__.__name__), "model results")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("radaroverlay", "[%s %s]" % ("1x1", obj.radaroverlay.__class__.__name__), "radar image for plot overlay")) + s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields")) + return s # }}} Index: /issm/trunk-jpl/src/m/classes/organizer.m =================================================================== --- /issm/trunk-jpl/src/m/classes/organizer.m (revision 25757) +++ /issm/trunk-jpl/src/m/classes/organizer.m (revision 25758) @@ -156,4 +156,26 @@ error(['Could not find ' path ]); end%}}} + function loaddatanoprefix(org,string),% {{{ + + %Get model path + if ~ischar(string), error('argument provided is not a string'); end + path=[org.repository '/' string]; + + %figure out if the data is there, otherwise, we have to use the default path supplied by user. + if exist(path,'file'), + path=path; + elseif exist([path '.mat'],'file'), + path=[path '.mat']; + else + error(['Could not find ' path ]); + end + if exist(path,'file') + evalin('caller',['load -mat ' path]); + return; + end + + %If we are here, the data has not been found. + error(['Could not find ' path ]); + end%}}} function bool=perform(org,varargin) % {{{ @@ -254,4 +276,30 @@ eval(['save(''' name '''' variables ',''-v7.3'');']); end%}}} + function savedatanoprefix(org,varargin) % {{{ + + %check + if (org.currentstep==0), error('Cannot save data because organizer (org) is empty! Make sure you did not skip any perform call'); end + if (org.currentstep>length(org.steps)), error('Cannot save data because organizer (org) is not up to date!'); end + + name=[org.repository '/' org.steps(org.currentstep).string ]; + disp(['saving data in: ' name]); + + %Skip if requested + if org.skipio, + disp(['WARNING: Skipping saving ' name]); + return; + end + + %check that md is a model + if (org.currentstep>length(org.steps)), error(['organizer error message: element with id ' num2str(org.currentstep) ' not found']); end + + %list of variable names: + variables=''; + for i=2:nargin, + variables=[variables ',' '''' inputname(i) '''']; + eval([inputname(i) '= varargin{' num2str(i-1) '};']); + end + eval(['save(''' name '''' variables ',''-v7.3'');']); + end%}}} end end Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25757) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25758) @@ -14,6 +14,6 @@ requested_outputs = {}; transitions = {}; - partitionice = []; - partitionhydro = []; + partitionice = []; + partitionhydro = []; end methods @@ -21,12 +21,10 @@ switch nargin case 0 - self=setdefaultparameters(self,earth); - case 1 - self=setdefaultparameters(self,varargin{:}); + self=setdefaultparameters(self); otherwise - error('solidearth constructor error message: zero or one argument only!'); + error('constructor not supported'); end end % }}} - function self = setdefaultparameters(self,planet) % {{{ + function self = setdefaultparameters(self) % {{{ %output default: @@ -35,5 +33,5 @@ %transitions should be a cell array of vectors: self.transitions={}; - + %no partitions requested for barystatic contribution: self.partitionice=[]; @@ -41,5 +39,5 @@ %earth radius - self.planetradius= planetradius(planet); + self.planetradius= planetradius('earth'); end % }}} Index: /issm/trunk-jpl/src/m/classes/solidearth.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.py (revision 25757) +++ /issm/trunk-jpl/src/m/classes/solidearth.py (revision 25758) @@ -1,3 +1,4 @@ import numpy as np + from checkfield import checkfield from fielddisplay import fielddisplay @@ -19,12 +20,14 @@ def __init__(self, *args): #{{{ - self.sealevel = np.nan - self.settings = solidearthsettings() - self.surfaceload = surfaceload() - self.lovenumbers = lovenumbers() - self.rotational = rotational() - self.planetradius = planetradius('earth') - self.requested_outputs = ['default'] - self.transitions = [] + self.sealevel = np.nan + self.settings = solidearthsettings() + self.surfaceload = surfaceload() + self.lovenumbers = lovenumbers() + self.rotational = rotational() + self.planetradius = planetradius('earth') + self.requested_outputs = [] + self.transitions = [] + self.partitionice = [] + self.partitionhydro = [] nargin = len(args) @@ -50,5 +53,16 @@ def setdefaultparameters(self): # {{{ - return self + # Default output + self.requested_outputs = ['default'] + + # Transitions should be a list + self.transitions = [] + + # No partitions requested for barystatic contribution + self.partitionice = [] + self.partitionhydro = [] + + # Earth radius + self.planetradius = planetradius('earth') #}}} @@ -58,4 +72,5 @@ md = checkfield(md, 'fieldname', 'solidearth.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) md = checkfield(md, 'fieldname', 'solidearth.requested_outputs', 'stringrow', 1) + self.settings.checkconsistency(md, solution, analyses) self.surfaceload.checkconsistency(md, solution, analyses) @@ -73,4 +88,19 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'planetradius', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'fieldname', 'transitions', 'format', 'MatArray') + + if len(self.partitionice): + npartice = np.max(self.partitionice) + 2 + else: + npartice = 0 + + if len(self.partitionhydro): + nparthydro = np.max(self.partitionhydro) + 2 + else: + nparthydro = 0 + + WriteData(fid, prefix, 'object', self, 'fieldname', 'partitionice', 'mattype', 1, 'format', 'DoubleMat'); + WriteData(fid, prefix, 'data', npartice, 'format', 'Integer', 'name', 'md.solidearth.npartice'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'partitionhydro', 'mattype', 1, 'format', 'DoubleMat'); + WriteData(fid, prefix, 'data', nparthydro,'format', 'Integer', 'name','md.solidearth.nparthydro'); self.settings.marshall(prefix, md, fid) Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25757) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25758) @@ -17,4 +17,5 @@ degacc = 0; %degree increment for resolution of Green tables horiz = 0; %compute horizontal deformation + glfraction = 1; %barystatic contribution full or fractional (default fractional) end methods @@ -48,4 +49,7 @@ %how many time steps we skip before we run solidearthsettings solver during transient self.runfrequency=1; + + %fractional contribution: + self.glfraction=1; %horizontal displacement? (not by default) @@ -63,9 +67,5 @@ md = checkfield(md,'fieldname','solidearth.settings.degacc','size',[1 1],'>=',1e-10); md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]); - - %checks on computational flags: - if self.elastic==1 & self.rigid==0, - error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set'); - end + md = checkfield(md,'fieldname','solidearth.settings.glfraction','values',[0 1]); %a coupler to a planet model is provided. @@ -96,4 +96,5 @@ fielddisplay(self,'rotation','earth rotational potential perturbation'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); + fielddisplay(self,'glfraction','contribute fractionally (default, 1) to barystatic sea level'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -109,4 +110,5 @@ WriteData(fid,prefix,'object',self,'fieldname','horiz','name','md.solidearth.settings.horiz','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','computesealevelchange','name','md.solidearth.settings.computesealevelchange','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','glfraction','name','md.solidearth.settings.glfraction','format','Integer'); end % }}} function savemodeljs(self,fid,modelname) % {{{ @@ -121,4 +123,5 @@ writejsdouble(fid,[modelname '.slr.settings.run_frequency'],self.run_frequency); writejsdouble(fid,[modelname '.slr.settings.degacc'],self.degacc); + writejsdouble(fid,[modelname '.slr.settings.glfraction'],self.glfraction); end % }}} function self = extrude(self,md) % {{{ Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25757) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25758) @@ -14,4 +14,17 @@ def __init__(self, *args): #{{{ + self.reltol = 0 + self.abstol = 0 + self.maxiter = 0 + self.rigid = 0 + self.elastic = 0 + self.rotation = 0 + self.ocean_area_scaling = 0 + self.runfrequency = 1 # How many time steps we skip before we run grd_core + self.computesealevelchange = 0 # Will grd_core compute sea level? + self.degacc = 0 # Degree increment for resolution of Green tables + self.horiz = 0 # Compute horizontal displacement? + self.glfraction = 1 # Barystatic contribution: full or fractional (default: fractional) + nargin = len(args) @@ -33,4 +46,5 @@ s += '{}\n'.format(fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation')) s += '{}\n'.format(fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green\'s functions')) + s += '{}\n'.format(fielddisplay(self, 'glfraction', 'contribute fractionally (default, 1) to barystatic sea level')) return s #}}} @@ -59,5 +73,7 @@ # Horizontal displacement? (not on by default) self.horiz = 0 - return self + + # Fractional contribution + self.glfraction = 1 #}}} @@ -71,4 +87,5 @@ md = checkfield(md, 'fieldname', 'solidearth.settings.degacc', 'size', [1], '>=', 1e-10) md = checkfield(md, 'fieldname', 'solidearth.settings.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'solidearth.settings.glfraction', 'values', [0, 1]) # A coupler to planet model is provided @@ -88,14 +105,15 @@ def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double') - WriteData(fid,prefix,'object', self,'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double') - WriteData(fid,prefix,'object', self,'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer') - WriteData(fid,prefix,'object', self,'fieldname', 'rigid', 'name', 'md.solidearth.settings.rigid', 'format', 'Boolean') - WriteData(fid,prefix,'object', self,'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean') - WriteData(fid,prefix,'object', self,'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean') - WriteData(fid,prefix,'object', self,'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean') - WriteData(fid,prefix,'object', self,'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer') - WriteData(fid,prefix,'object', self,'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double') - WriteData(fid,prefix,'object', self,'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer') - WriteData(fid,prefix,'object', self,'fieldname', 'computesealevelchange', 'name', 'md.solidearth.settings.computesealevelchange', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'rigid', 'name', 'md.solidearth.settings.rigid', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'computesealevelchange', 'name', 'md.solidearth.settings.computesealevelchange', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname','glfraction', 'name', 'md.solidearth.settings.glfraction', 'format', 'Integer') #}}} Index: /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m (revision 25757) +++ /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m (revision 25758) @@ -183,4 +183,5 @@ %Go through O2 regions: for i=subsetregions, + %for i=33, string=self.boxes(i).RGI_CODE; disp(['progressing with region ' num2str(i) ' ' string]); @@ -208,5 +209,5 @@ case 19, radius=60; case 32, radius=60; - case 33, radius=5; + case 33, radius=10; case 41, radius=75; case 42, radius=45; Index: /issm/trunk-jpl/src/m/contrib/larour/ismip6.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/ismip6.m (revision 25757) +++ /issm/trunk-jpl/src/m/contrib/larour/ismip6.m (revision 25758) @@ -13,6 +13,12 @@ directories = {}; %directories where the files are experiments = {}; %names of experiments + base = {}; %placeholder for base + surface = {}; %placeholder for surface thickness = {}; %placeholder for thicknesses deltathickness = {}; %placeholder for delta thicknesses + deltathicknessvaf = {}; %placeholder for delta thicknesses above floatation + deltathicknesshal = {}; %placeholder for delta thicknesses halosteric origins + deltathicknessbar = {}; %placeholder for delta thicknesses halosteric origins + thicknesscorrection={}; icemask = {}; %placeholder for ice masks oceanmask = {}; %placeholder for ocean masks @@ -20,4 +26,5 @@ timestart = {}; %placeholder for times calendar = {}; %placeholder for times + di = {}; %ice densities end methods @@ -146,4 +153,7 @@ end + if ~exist('timestart','var'), timestart=2015; end + if ~exist('calendar','var'), calendar=0; end + end % }}} function info=readinfo(self,experiment,field) % {{{ @@ -203,5 +213,10 @@ %map onto mesh: correct only for thicknesses - hg=ismip2mesh_correction.*(ismip2mesh*ht) ; + if strcmpi(field,'lithk') | strcmpi(field,'orog') | strcmpi(field,'base'), + hg=ismip2mesh_correction.*(ismip2mesh*ht) ; + %hg=ismip2mesh*ht ; + else + hg=ismip2mesh*ht ; + end %keep field: @@ -209,4 +224,12 @@ pos=find(isnan(hg)); hg(pos)=0; self.thickness{i}=hg; + end + if strcmpi(field,'orog'), + pos=find(isnan(hg)); hg(pos)=0; + self.surface{i}=hg; + end + if strcmpi(field,'base'), + pos=find(isnan(hg)); hg(pos)=0; + self.base{i}=hg; end if strcmpi(field,'sftgif'), @@ -220,11 +243,12 @@ end if strcmpi(field,'sftgrf'), - hge=-ones(md.mesh.numberofvertices,size(hg,2)); + hgv=-ones(md.mesh.numberofvertices,size(hg,2)); for j=1:size(hg,2), hgj=hg(:,j); - pos=find(hgj>0); - hge(md.mesh.elements(pos,:),j)=1; - end - self.oceanmask{i}=hge; + pos=find(hgj>.99); %we want fully grounded + %pos=find(hgj>0); %we want slightly grounded + hgv(md.mesh.elements(pos,:),j)=1; + end + self.oceanmask{i}=hgv; end Index: /issm/trunk-jpl/src/m/contrib/larour/mme_autotime_correlation_matrix.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/mme_autotime_correlation_matrix.m (revision 25757) +++ /issm/trunk-jpl/src/m/contrib/larour/mme_autotime_correlation_matrix.m (revision 25758) @@ -1,53 +1,20 @@ -function matrix=mme_time_correlation_matrix(mme,varargin); +function matrix=mme_autotime_correlation_matrix(mme,type) - if nargin==2, - type=varargin{1}; - elseif nargin==3, - mme2=varargin{1}; - type=varargin{2}; - else - error('mme_time_correlation_matrix usage error: 2 or 3 arguments only allowed!'); + + %Out of a multi model ensemble (nsamples x nsteps) of runs, build + %a temporal correlation matrix (of size nsteps x nsteps) + + nsamples=size(mme,1); + nsteps=size(mme,2); + + %initialize with 1 in the diagonal: + matrix=eye(nsteps,nsteps); + + %go through time steps, and fill up the top part. + for i=1:nsteps, + for j=i+1:nsteps, + matrix(i,j)=corr(mme(:,i),mme(:,j),'Type',type); + matrix(j,i)=matrix(i,j); + end end - if nargin==2, - - %Out of a multi model ensemble (nsamples x nsteps) of runs, build - %a temporal correlation matrix (of size nsteps x nsteps) - - nsamples=size(mme,1); - nsteps=size(mme,2); - - %initialize with 1 in the diagonal: - matrix=eye(nsteps,nsteps); - - %go through time steps, and fill up the top part. - for i=1:nsteps, - for j=i+1:nsteps, - matrix(i,j)=corr(mme(:,i),mme(:,j),'Type',type); - matrix(j,i)=matrix(i,j); - end - end - else - - %Same kind of computations, except it's not autocorrelation: - nsamples=size(mme,1); nsamples2=size(mme2,1); - nsteps=size(mme,2); nsteps2=size(mme2,2); - - if nsteps2~=nsteps, - error('number of time steps from both sample matrices should be identical!'); - end - if nsamples2~=nsamples, - error('number of samples from both sample matrices should be identical!'); - end - - %initialize with 1 in the diagonal: - matrix=zeros(nsteps,nsteps); - - %go through time steps, and fill up the top part. - for i=1:nsteps, - for j=i:nsteps, - matrix(i,j)=corr(mme(:,i),mme2(:,j),'Type',type); - matrix(j,i)=matrix(i,j); - end - end - end Index: /issm/trunk-jpl/src/m/interp/averaging.m =================================================================== --- /issm/trunk-jpl/src/m/interp/averaging.m (revision 25757) +++ /issm/trunk-jpl/src/m/interp/averaging.m (revision 25758) @@ -64,5 +64,9 @@ else rep=3; - areas=GetAreas(index,md.mesh.x2d,md.mesh.y2d); + if isa(md.mesh,'mesh3dsurface'), + areas=GetAreas3DTria(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z); + else + areas=GetAreas(index,md.mesh.x,md.mesh.y); + end end summation=1/rep*ones(rep,1); Index: /issm/trunk-jpl/src/m/mesh/FixMesh.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/FixMesh.m (revision 25757) +++ /issm/trunk-jpl/src/m/mesh/FixMesh.m (revision 25758) @@ -1,3 +1,3 @@ -function [index2 x2 y2 value2]=FixMesh(index,x,y,value) +function [index2 x2 y2 value2 newpos]=FixMesh(index,x,y,value) % FIXMESH - FixMesh fix mesh with broken triangles, orphan vertices, etc ... % @@ -15,4 +15,5 @@ y2=y; value2=value; +newpos=1:length(x); %First, look for orphan vertices, and take them out. @@ -29,4 +30,5 @@ y2(orphan)=[]; value2(orphan)=[]; + newpos(orphan)=[]; %now, the index: Index: /issm/trunk-jpl/src/m/plot/processdatalatlong.m =================================================================== --- /issm/trunk-jpl/src/m/plot/processdatalatlong.m (revision 25757) +++ /issm/trunk-jpl/src/m/plot/processdatalatlong.m (revision 25758) @@ -22,8 +22,6 @@ %interpolate data: extradata=griddata(x0,y0,data,xextra,yextra,'nearest'); - data=[data; extradata]; elseif length(data)==length(md.mesh.elements), - error('processdatalatlong error message: coord ''latlong'' case not covered for element data '); datatype=1; end Index: /issm/trunk-jpl/src/m/solve/solve.m =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.m (revision 25757) +++ /issm/trunk-jpl/src/m/solve/solve.m (revision 25758) @@ -4,29 +4,31 @@ % Usage: % md=solve(md,solutionstring,varargin) -% where varargin is a lit of paired arguments of string OR enums % -% solution types available comprise: -% - 'Stressbalance' or 'sb' -% - 'Masstransport' or 'mt' -% - 'Thermal' or 'th' -% - 'Steadystate' or 'ss' -% - 'Transient' or 'tr' -% - 'Balancethickness' or 'mc' -% - 'Balancevelocity' or 'bv' -% - 'BedSlope' or 'bsl' -% - 'SurfaceSlope' or 'ssl' -% - 'Hydrology' or 'hy' -% - 'DamageEvolution' or 'da' -% - 'Gia' or 'gia' -% - 'Love' or 'lv' -% - 'Esa' or 'esa' -% - 'Sealevelrise' or 'slr' +% where varargin is a lit of paired arguments of string OR enums % -% extra options: -% - loadonly : does not solve. only load results -% - runtimename : true or false (default is true), makes name unique -% - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on consistency of model -% - restart: 'directory name (relative to the execution directory) where the restart file is located. -% - outbinread : if 0, download the outbin but do not process is (md.results is not updated) +% Solution types available comprise: +% - 'Stressbalance' or 'sb' +% - 'Masstransport' or 'mt' +% - 'Thermal' or 'th' +% - 'Steadystate' or 'ss' +% - 'Transient' or 'tr' +% - 'Balancethickness' or 'mc' +% - 'Balancevelocity' or 'bv' +% - 'BedSlope' or 'bsl' +% - 'SurfaceSlope' or 'ssl' +% - 'Hydrology' or 'hy' +% - 'DamageEvolution' or 'da' +% - 'Gia' or 'gia' +% - 'Love' or 'lv' +% - 'Esa' or 'esa' +% - 'Sealevelrise' or 'slr' +% +% Extra options: +% - loadonly : does not solve. only load results +% - runtimename : true or false (default is true), makes name unique +% - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on +% consistency of model +% - restart : 'directory name (relative to the execution directory) +% where the restart file is located % % Examples: @@ -111,5 +113,5 @@ end -%if running qmu analysis, some preprocessing of dakota files using models fields needs to be carried out. +%if running QMU analysis, some preprocessing of Dakota files using models fields needs to be carried out. if md.qmu.isdakota, md=preqmu(md,options); Index: /issm/trunk-jpl/src/m/solve/solve.py =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.py (revision 25757) +++ /issm/trunk-jpl/src/m/solve/solve.py (revision 25758) @@ -1,3 +1,3 @@ -import datetime +from datetime import datetime import os import shutil @@ -16,29 +16,31 @@ Usage: md = solve(md, solutionstring, varargin) - where varargin is a list of paired arguments of string OR enums + + where varargin is a list of paired arguments of string OR enums - solution types available comprise: - - 'Stressbalance' or 'sb' - - 'Masstransport' or 'mt' - - 'Thermal' or 'th' - - 'Steadystate' or 'ss' - - 'Transient' or 'tr' - - 'Balancethickness' or 'mc' - - 'Balancevelocity' or 'bv' - - 'BedSlope' or 'bsl' - - 'SurfaceSlope' or 'ssl' - - 'Hydrology' or 'hy' - - 'DamageEvolution' or 'da' - - 'Gia' or 'gia' - - 'Esa' or 'esa' - - 'Sealevelrise' or 'slr' - - 'Love' or 'lv' + Solution types available comprise: + - 'Stressbalance' or 'sb' + - 'Masstransport' or 'mt' + - 'Thermal' or 'th' + - 'Steadystate' or 'ss' + - 'Transient' or 'tr' + - 'Balancethickness' or 'mc' + - 'Balancevelocity' or 'bv' + - 'BedSlope' or 'bsl' + - 'SurfaceSlope' or 'ssl' + - 'Hydrology' or 'hy' + - 'DamageEvolution' or 'da' + - 'Gia' or 'gia' + - 'Esa' or 'esa' + - 'Sealevelrise' or 'slr' + - 'Love' or 'lv' - extra options: - - loadonly : does not solve. only load results - - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures - checks on consistency of model - - restart : 'directory name (relative to the execution - directory) where the restart file is located + Extra options: + - loadonly : does not solve. only load results + - runtimename : true or false (default is true), makes name unique + - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on + consistency of model + - restart : directory name (relative to the execution directory) + where the restart file is located Examples: @@ -47,5 +49,5 @@ """ - #recover and process solve options + # Recover and process solve options if solutionstring.lower() == 'sb' or solutionstring.lower() == 'stressbalance': solutionstring = 'StressbalanceSolution' @@ -82,5 +84,5 @@ options = pairoptions('solutionstring', solutionstring, *args) - #recover some fields + # Recover some fields md.private.solution = solutionstring cluster = md.cluster @@ -90,5 +92,5 @@ batch = 0 - #check model consistency + # Check model consistency if options.getfieldvalue('checkconsistency', 'yes') == 'yes': if md.verbose.solution: @@ -96,8 +98,8 @@ ismodelselfconsistent(md) - #First, build a runtime name that is unique + # First, build a runtime name that is unique restart = options.getfieldvalue('restart', '') if restart == 1: - pass #do nothing + pass # do nothing else: if restart: @@ -105,25 +107,25 @@ else: if options.getfieldvalue('runtimename', True): - c = datetime.datetime.now() + c = datetime.now() md.private.runtimename = "%s-%02i-%02i-%04i-%02i-%02i-%02i-%i" % (md.miscellaneous.name, c.month, c.day, c.year, c.hour, c.minute, c.second, os.getpid()) else: md.private.runtimename = md.miscellaneous.name - #if running qmu analysis, some preprocessing of dakota files using models - #fields needs to be carried out. + # If running QMU analysis, some preprocessing of Dakota files using models + # fields needs to be carried out if md.qmu.isdakota: md = preqmu(md, options) - #Do we load results only? + # Do we load results only? if options.getfieldvalue('loadonly', False): md = loadresultsfromcluster(md) return md - #Write all input files - marshall(md) # bin file - md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits') # toolkits file - cluster.BuildQueueScript(md.private.runtimename, md.miscellaneous.name, md.private.solution, md.settings.io_gather, md.debug.valgrind, md.debug.gprof, md.qmu.isdakota, md.transient.isoceancoupling) # queue file + # Write all input files + marshall(md) # bin file + md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits') # toolkits file + cluster.BuildQueueScript(md.private.runtimename, md.miscellaneous.name, md.private.solution, md.settings.io_gather, md.debug.valgrind, md.debug.gprof, md.qmu.isdakota, md.transient.isoceancoupling) # queue file - #Stop here if batch mode + # Stop here if batch mode if options.getfieldvalue('batch', 'no') == 'yes': print('batch mode requested: not launching job interactively') @@ -131,5 +133,5 @@ return md - #Upload all required files: + # Upload all required files: modelname = md.miscellaneous.name filelist = [modelname + '.bin ', modelname + '.toolkits ', modelname + '.queue '] @@ -140,16 +142,16 @@ cluster.UploadQueueJob(md.miscellaneous.name, md.private.runtimename, filelist) - #Launch job + # Launch job cluster.LaunchQueueJob(md.miscellaneous.name, md.private.runtimename, filelist, restart, batch) - #wait on lock + # Wait on lock if md.settings.waitonlock > 0: - #we wait for the done file islock = waitonlock(md) - if islock == 0: #no results to be loaded + if islock == 0: # no results to be loaded print('The results must be loaded manually with md = loadresultsfromcluster(md).') - else: #load results + else: # load results if md.verbose.solution: print('loading results from cluster') md = loadresultsfromcluster(md) + return md Index: /issm/trunk-jpl/src/m/units/sletogt.m =================================================================== --- /issm/trunk-jpl/src/m/units/sletogt.m (revision 25757) +++ /issm/trunk-jpl/src/m/units/sletogt.m (revision 25758) @@ -1,3 +1,4 @@ function conversionfactor=sletogt() - - conversionfactor=361.9; %361.9 Gigatons to 1 mm sea-level equivalent. + + rho_water=1023; + conversionfactor=rho_water/1000*361.9; %361.9 Gigatons to 1 mm sea-level equivalent.