Index: /issm/trunk-jpl/externalpackages/gmsh/install-4.sh =================================================================== --- /issm/trunk-jpl/externalpackages/gmsh/install-4.sh (revision 25766) +++ /issm/trunk-jpl/externalpackages/gmsh/install-4.sh (revision 25767) @@ -6,4 +6,6 @@ # VER="4.5.6" + +PETSC_ROOT="${ISSM_DIR}/externalpackages/petsc/install" # Cleanup @@ -41,5 +43,6 @@ -DENABLE_MPI=1 \ -DENABLE_OCC=0 \ - -DENABLE_TOUCHBAR=0 + -DENABLE_TOUCHBAR=0 \ + -DMETIS_ROOT="${PETSC_ROOT}" # Compile and install Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25766) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25767) @@ -5446,6 +5446,6 @@ /*Computational flags:*/ - bool computerigid = true; - bool computeelastic = true; + bool computerigid = false; + bool computeelastic = false; int horiz; @@ -5469,4 +5469,8 @@ parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); + /*allocate indices:*/ + indices=xNew(gsize); + G=xNewZeroInit(gsize); + if(computeelastic){ parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); @@ -5478,10 +5482,5 @@ 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){ + GU=xNewZeroInit(gsize); if(horiz){ @@ -5549,6 +5548,6 @@ /*Add in inputs:*/ - this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); - this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); + this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); + this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); if(computeelastic){ this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); @@ -5558,4 +5557,6 @@ } } + this->inputs->SetDoubleInput(AreaEnum,this->lid,area); + this->AddInput(SealevelAreaEnum,&area,P0Enum); /*Free allocations:*/ @@ -5594,5 +5595,5 @@ bool notfullygrounded=false; bool scaleoceanarea= false; - bool computeelastic= true; + bool computerigid= false; int glfraction=1; @@ -5640,10 +5641,10 @@ rho_ice=FindParam(MaterialsRhoIceEnum); rho_water=FindParam(MaterialsRhoSeawaterEnum); - this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); + this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); /*retrieve precomputed G:*/ - if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); + if(computerigid)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ @@ -5699,5 +5700,5 @@ _assert_(!xIsNan(bslrice)); - if(computeelastic){ + if(computerigid){ /*convert from m to kg/m^2:*/ I=I*rho_ice*phi; @@ -5725,5 +5726,5 @@ bool notfullygrounded=false; bool scaleoceanarea= false; - bool computeelastic= true; + bool computeelastic= false; /*elastic green function:*/ @@ -5799,5 +5800,5 @@ IssmDouble BP; //change in bottom pressure (Farrel and Clarke, Equ. 4) IssmDouble constant; - bool computeelastic= true; + bool computeelastic= false; /*elastic green function:*/ @@ -5901,5 +5902,5 @@ /*computational flags:*/ - bool computeelastic= true; + bool computeelastic= false; /*early return if we are not on the ocean or on an ice cap:*/ Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp (revision 25766) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp (revision 25767) @@ -12,5 +12,5 @@ #include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h" #include "../InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h" - + void InputUpdateFromDakotax(FemModel* femmodel,double* variables,char* *variables_descriptors,int numdakotavariables){ /*{{{*/ @@ -34,27 +34,27 @@ /*retrieve parameters: */ - femmodel->parameters->FindParam(&variable_partitions,&variable_partitions_num,NULL,NULL,QmuVariablePartitionsEnum); - femmodel->parameters->FindParam(&variable_partitions_npart,NULL,NULL,QmuVariablePartitionsNpartEnum); - femmodel->parameters->FindParam(&variable_partitions_nt,NULL,NULL,QmuVariablePartitionsNtEnum); - - - /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and - * for each descriptor, take the variable value and plug it into the inputs (more or less :)): - * We also start with distributed and standard values , as they tend to be used to pluck data from a multi-modle ensemble (mme) - * which can then be scaled. Doing the scaling first would be impractical, as the entire mme would have to be scaled, + femmodel->parameters->FindParam(&variable_partitions,&variable_partitions_num,NULL,NULL,QmuVariablePartitionsEnum); + femmodel->parameters->FindParam(&variable_partitions_npart,NULL,NULL,QmuVariablePartitionsNpartEnum); + femmodel->parameters->FindParam(&variable_partitions_nt,NULL,NULL,QmuVariablePartitionsNtEnum); + + + /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and + * for each descriptor, take the variable value and plug it into the inputs (more or less :)): + * We also start with distributed and standard values , as they tend to be used to pluck data from a multi-modle ensemble (mme) + * which can then be scaled. Doing the scaling first would be impractical, as the entire mme would have to be scaled, * which is a waste of time:*/ variablecount=0; - for(i=0;i DragCoefficient). + /*Variable is distributed. Determine root name of variable (ex: distributed_DragCoefficient_1 -> DragCoefficient). * Allocate distributed_values and fill the distributed_values with the next npart variables: */ @@ -86,11 +86,11 @@ //for (int j=0;j(parameter); @@ -104,16 +104,16 @@ variablecount++; } - + variablecount=0; /*now deal with scaled variabes:*/ - for(i=0;i DragCoefficient). Allocate distributed_values and fill the + /* Variable is scaled, determine its root name (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the * distributed_values with the next npart variables coming from Dakota: */ memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char)); @@ -136,5 +136,5 @@ /*increment i to skip the distributed values just collected: */ i+=npart*nt-1; //careful, the for loop will add 1. - + /*Free allocations: */ xDelete(parameter); @@ -202,44 +202,33 @@ transientinput2 = femmodel->inputs->GetTransientInput(DummyEnum); - for (int p=npart;p>=0;p--){ - int pp=p; - if (p==npart)pp=-1; /*so, the logic is, we want to do the -1 partition first, then - go from npart-1 to 0 in reverse order:*/ - - for (int i=0;ielements->Size();i++){ - int element_partition; - - Tria* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - - element_partition= (int)variable_partition[element->Sid()]; - if(element_partition!=pp)continue; - - if(element_partition==-1)id=0; //grab background field - else id=distributed_values[element_partition]-1; //grab partition field - - /*recover the right field from the mme: */ - transientinput = datasetinput->GetTransientInputByOffset(id); _assert_(transientinput); - - /*copy values from the transientinput to the final transientinput2: */ - for (int j=0;jGetTriaInput(j); - element->InputServe(tria_input); - if(interpolationenum==P0Enum){ - value=tria_input->element_values[0]; - transientinput2->AddTriaTimeInput( j,1,&(element->lid),&value,P0Enum); - } - else if(interpolationenum==P1Enum){ - - /*Get values and lid list*/ - const int numvertices = element->GetNumberOfVertices(); - int *vertexlids = xNew(numvertices); - int *vertexsids = xNew(numvertices); - - /*Recover vertices ids needed to initialize inputs*/ - element->GetVerticesLidList(&vertexlids[0]); - element->GetVerticesSidList(&vertexsids[0]); - values=tria_input->element_values; - transientinput2->AddTriaTimeInput( j,numvertices,vertexlids,values,P1Enum); - } + for(Object* & object : femmodel->elements->objects){ + Tria* element=xDynamicCast(object); + + if((int)variable_partition[element->Sid()]==-1)id=0; //grab background field + else id=distributed_values[(int)variable_partition[element->Sid()]]-1; //grab partition field + + /*recover the right field from the mme: */ + transientinput = datasetinput->GetTransientInputByOffset(id); _assert_(transientinput); + + /*copy values from the transientinput to the final transientinput2: */ + for (int j=0;jGetTriaInput(j); + element->InputServe(tria_input); + if(interpolationenum==P0Enum){ + value=tria_input->element_values[0]; + transientinput2->AddTriaTimeInput( j,1,&(element->lid),&value,P0Enum); + } + else if(interpolationenum==P1Enum){ + + /*Get values and lid list*/ + const int numvertices = element->GetNumberOfVertices(); + int *vertexlids = xNew(numvertices); + int *vertexsids = xNew(numvertices); + + /*Recover vertices ids needed to initialize inputs*/ + element->GetVerticesLidList(&vertexlids[0]); + element->GetVerticesSidList(&vertexsids[0]); + values=tria_input->element_values; + transientinput2->AddTriaTimeInput( j,numvertices,vertexlids,values,P1Enum); } } @@ -258,10 +247,10 @@ /*Go through elements, copy input name to dummy, and scale it using the distributed_values and the partition vector:*/ - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); element->InputScaleFromDakota(distributed_values,partition,npart,nt,name); } - /*We created a dummy input, which was a scaled copy of the name input. Now wipe + /*We created a dummy input, which was a scaled copy of the name input. Now wipe * out the name input with the new input:*/ femmodel->inputs->ChangeEnum(DummyEnum,name); Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 25766) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 25767) @@ -195,12 +195,59 @@ switch nargin case 0 - md=setdefaultparameters(md); - case 1 - error('model constructor not supported yet'); - + md=setdefaultparameters(md,'earth'); otherwise - error('model constructor error message: 0 of 1 argument only in input.'); - end - + options=pairoptions(varargin{:}); + planet=getfieldvalue(options,'planet','earth'); + md=setdefaultparameters(md,planet); + end + + end + %}}} + function md = setdefaultparameters(md,planet) % {{{ + + %initialize subclasses + md.mesh = mesh2d(); + md.mask = mask(); + md.constants = constants(); + md.geometry = geometry(); + md.initialization = initialization(); + md.smb = SMBforcing(); + md.basalforcings = basalforcings(); + md.friction = friction(); + md.rifts = rifts(); + md.solidearth = solidearth(planet); + md.dsl = dsl(); + md.timestepping = timestepping(); + md.groundingline = groundingline(); + md.materials = matice(); + md.damage = damage(); + md.flowequation = flowequation(); + md.debug = debug(); + md.verbose = verbose(); + md.settings = issmsettings(); + md.toolkits = toolkits(); + md.cluster = generic(); + md.balancethickness = balancethickness(); + md.stressbalance = stressbalance(); + md.hydrology = hydrologyshreve(); + md.masstransport = masstransport(); + md.thermal = thermal(); + md.steadystate = steadystate(); + md.transient = transient(); + md.levelset = levelset(); + md.calving = calving(); + md.frontalforcings = frontalforcings(); + md.gia = giamme(); + md.love = fourierlove(); + md.esa = esa(); + md.autodiff = autodiff(); + md.inversion = inversion(); + md.qmu = qmu(); + md.amr = amr(); + md.radaroverlay = radaroverlay(); + md.results = struct(); + md.outputdefinition = outputdefinition(); + md.miscellaneous = miscellaneous(); + md.private = private(); end %}}} @@ -1288,52 +1335,4 @@ end end% }}} - function md = setdefaultparameters(md) % {{{ - - %initialize subclasses - md.mesh = mesh2d(); - md.mask = mask(); - md.constants = constants(); - md.geometry = geometry(); - md.initialization = initialization(); - md.smb = SMBforcing(); - md.basalforcings = basalforcings(); - md.friction = friction(); - md.rifts = rifts(); - md.solidearth = solidearth(); - md.dsl = dsl(); - md.timestepping = timestepping(); - md.groundingline = groundingline(); - md.materials = matice(); - md.damage = damage(); - md.flowequation = flowequation(); - md.debug = debug(); - md.verbose = verbose(); - md.settings = issmsettings(); - md.toolkits = toolkits(); - md.cluster = generic(); - md.balancethickness = balancethickness(); - md.stressbalance = stressbalance(); - md.hydrology = hydrologyshreve(); - md.masstransport = masstransport(); - md.thermal = thermal(); - md.steadystate = steadystate(); - md.transient = transient(); - md.levelset = levelset(); - md.calving = calving(); - md.frontalforcings = frontalforcings(); - md.gia = giamme(); - md.love = fourierlove(); - md.esa = esa(); - md.autodiff = autodiff(); - md.inversion = inversion(); - md.qmu = qmu(); - md.amr = amr(); - md.radaroverlay = radaroverlay(); - md.results = struct(); - md.outputdefinition = outputdefinition(); - md.miscellaneous = miscellaneous(); - md.private = private(); - end - %}}} function md = tetras(md,varargin) % {{{ %TETRAS - split 3d prismatic mesh into 3 tetrahedrons Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 25766) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 25767) @@ -3,4 +3,5 @@ import copy import sys +from pairoptions import pairoptions from mesh2d import mesh2d from mesh3dprisms import mesh3dprisms @@ -79,49 +80,59 @@ class model(object): #properties - def __init__(self): #{{{ - self.mesh = mesh2d() - self.mask = mask() - self.constants = constants() - self.geometry = geometry() - self.initialization = initialization() - self.smb = SMBforcing() - self.basalforcings = basalforcings() - self.friction = friction() - self.rifts = rifts() - self.solidearth = solidearth() - self.dsl = dsl() - self.timestepping = timestepping() - self.groundingline = groundingline() - self.materials = matice() - self.damage = damage() - self.flowequation = flowequation() - self.debug = debug() - self.verbose = verbose() - self.settings = issmsettings() - self.toolkits = toolkits() - self.cluster = generic() - self.balancethickness = balancethickness() - self.stressbalance = stressbalance() - self.hydrology = hydrologyshreve() - self.masstransport = masstransport() - self.thermal = thermal() - self.steadystate = steadystate() - self.transient = transient() - self.levelset = levelset() - self.calving = calving() - self.frontalforcings = frontalforcings() - self.gia = giamme() - self.love = fourierlove() - self.esa = esa() - self.autodiff = autodiff() - self.inversion = inversion() - self.qmu = qmu() - self.amr = amr() - self.radaroverlay = radaroverlay() - self.results = results() - self.outputdefinition = outputdefinition() - self.miscellaneous = miscellaneous() - self.private = private() - #}}} + def __init__(self, *args): #{{{ + self.mesh = None + self.mask = None + self.constants = None + self.geometry = None + self.initialization = None + self.smb = None + self.basalforcings = None + self.friction = None + self.rifts = None + self.solidearth = None + self.dsl = None + self.timestepping = None + self.groundingline = None + self.materials = None + self.damage = None + self.flowequation = None + self.debug = None + self.verbose = None + self.settings = None + self.toolkits = None + self.cluster = None + self.balancethickness = None + self.stressbalance = None + self.hydrology = None + self.masstransport = None + self.thermal = None + self.steadystate = None + self.transient = None + self.levelset = None + self.calving = None + self.frontalforcings = None + self.gia = None + self.love = None + self.esa = None + self.autodiff = None + self.inversion = None + self.qmu = None + self.amr = None + self.radaroverlay = None + self.results = None + self.outputdefinition = None + self.miscellaneous = None + self.private = None + + nargs = len(args) + + if nargs == 0: + self.setdefaultparameters('earth') + else: + self.setdefaultparameters(args[0]) + options = pairoptions(*args) + planet = options.getfieldvalue('planet', 'earth') + self.setdefaultparameters(planet) +#}}} def properties(self): # {{{ @@ -224,4 +235,50 @@ # }}} + def setdefaultparameters(self, planet): #{{{ + self.mesh = mesh2d() + self.mask = mask() + self.constants = constants() + self.geometry = geometry() + self.initialization = initialization() + self.smb = SMBforcing() + self.basalforcings = basalforcings() + self.friction = friction() + self.rifts = rifts() + self.solidearth = solidearth(planet) + self.dsl = dsl() + self.timestepping = timestepping() + self.groundingline = groundingline() + self.materials = matice() + self.damage = damage() + self.flowequation = flowequation() + self.debug = debug() + self.verbose = verbose() + self.settings = issmsettings() + self.toolkits = toolkits() + self.cluster = generic() + self.balancethickness = balancethickness() + self.stressbalance = stressbalance() + self.hydrology = hydrologyshreve() + self.masstransport = masstransport() + self.thermal = thermal() + self.steadystate = steadystate() + self.transient = transient() + self.levelset = levelset() + self.calving = calving() + self.frontalforcings = frontalforcings() + self.gia = giamme() + self.love = fourierlove() + self.esa = esa() + self.autodiff = autodiff() + self.inversion = inversion() + self.qmu = qmu() + self.amr = amr() + self.radaroverlay = radaroverlay() + self.results = results() + self.outputdefinition = outputdefinition() + self.miscellaneous = miscellaneous() + self.private = private() + #}}} + def checkmessage(self, string): # {{{ print("model not consistent: {}".format(string)) Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25766) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25767) @@ -6,5 +6,5 @@ classdef solidearth properties (SetAccess=public) - initialsealevel = NaN; + initialsealevel = NaN; settings = solidearthsettings(); external = []; @@ -22,26 +22,28 @@ switch nargin case 0 - self=setdefaultparameters(self); + self=setdefaultparameters(self,'earth'); + case 1 + self=setdefaultparameters(self,varargin{:}); otherwise - error('constructor not supported'); + error('solidearth constructor error message: zero or one argument only!'); end end % }}} - function self = setdefaultparameters(self) % {{{ - - %output default: - self.requested_outputs={'default'}; + function self = setdefaultparameters(self,planet) % {{{ - %transitions should be a cell array of vectors: - self.transitions={}; - - %no partitions requested for barystatic contribution: - self.partitionice=[]; - self.partitionhydro=[]; + %output default: + self.requested_outputs={'default'}; - %no external solutions by defalt: - self.external=[]; + %transitions should be a cell array of vectors: + self.transitions={}; + + %no partitions requested for barystatic contribution: + self.partitionice=[]; + self.partitionhydro=[]; - %earth radius - self.planetradius= planetradius('earth'); + %no external solutions by defalt: + self.external=[]; + + %earth radius + self.planetradius= planetradius(planet); end % }}} Index: /issm/trunk-jpl/src/m/classes/solidearth.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.py (revision 25766) +++ /issm/trunk-jpl/src/m/classes/solidearth.py (revision 25767) @@ -32,10 +32,12 @@ self.partitionhydro = [] - nargin = len(args) + nargs = len(args) - if nargin == 0: - self.setdefaultparameters() + if nargs == 0: + self.setdefaultparameters('earth') + elif nargs == 1: + self.setdefaultparameters(args[0]) else: - raise Exception('constructor not supported') + raise Exception('solidearth constructor error message: zero or one argument only!') #}}} @@ -59,5 +61,5 @@ #}}} - def setdefaultparameters(self): # {{{ + def setdefaultparameters(self, planet): # {{{ # Default output self.requested_outputs = ['default'] @@ -74,5 +76,5 @@ # Earth radius - self.planetradius = planetradius('earth') + self.planetradius = planetradius(planet) #}}} Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25766) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25767) @@ -71,4 +71,9 @@ md = checkfield(md,'fieldname','solidearth.settings.glfraction','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 + %a coupler to a planet model is provided. if self.isgrd, Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25766) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25767) @@ -92,4 +92,8 @@ md = checkfield(md, 'fieldname', 'solidearth.settings.glfraction', 'values', [0, 1]) + # Checks on computational flags + if self.elastic and not self.rigid: + raise Exception('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set') + # A coupler to planet model is provided if self.isgrd: Index: /issm/trunk-jpl/src/m/geometry/planetradius.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/planetradius.py (revision 25766) +++ /issm/trunk-jpl/src/m/geometry/planetradius.py (revision 25767) @@ -11,5 +11,7 @@ if planet == 'earth': - radius = 6.371012 * pow(10, 6) + radius = 6.371012e6 + elif planet == 'europa': + radius = 1.5008e6 else: raise TypeError("planet type %s not supported yet!" % planet)