Changeset 25767


Ignore:
Timestamp:
11/17/20 00:53:08 (4 years ago)
Author:
jdquinn
Message:

BUG: Fix for SE with rigid (corrected issue that I introduced earlier); various

Location:
issm/trunk-jpl
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/externalpackages/gmsh/install-4.sh

    r25746 r25767  
    66#
    77VER="4.5.6"
     8
     9PETSC_ROOT="${ISSM_DIR}/externalpackages/petsc/install"
    810
    911# Cleanup
    … …  
    4143        -DENABLE_MPI=1 \
    4244        -DENABLE_OCC=0 \
    43         -DENABLE_TOUCHBAR=0
     45        -DENABLE_TOUCHBAR=0 \
     46        -DMETIS_ROOT="${PETSC_ROOT}"
    4447
    4548# Compile and install

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

    r25763 r25767  
    54465446
    54475447        /*Computational flags:*/
    5448         bool computerigid = true;
    5449         bool computeelastic = true;
     5448        bool computerigid = false;
     5449        bool computeelastic = false;
    54505450        int  horiz;
    54515451
    … …  
    54695469        parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
    54705470
     5471        /*allocate indices:*/
     5472        indices=xNew<IssmDouble>(gsize);
     5473        G=xNewZeroInit<IssmDouble>(gsize);
     5474
    54715475        if(computeelastic){
    54725476                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
    … …  
    54785482                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
    54795483                parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
    5480         }
    5481 
    5482         /*allocate indices:*/
    5483         indices=xNew<IssmDouble>(gsize);
    5484         G=xNewZeroInit<IssmDouble>(gsize);
    5485         if(computeelastic){
     5484
    54865485                GU=xNewZeroInit<IssmDouble>(gsize);
    54875486                if(horiz){
    … …  
    55495548
    55505549        /*Add in inputs:*/
    5551     this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
    5552     this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
     5550        this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
     5551        this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
    55535552        if(computeelastic){
    55545553                this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
    … …  
    55585557                }
    55595558        }
     5559        this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
     5560        this->AddInput(SealevelAreaEnum,&area,P0Enum);
    55605561
    55615562        /*Free allocations:*/
    … …  
    55945595        bool notfullygrounded=false;
    55955596        bool scaleoceanarea= false;
    5596         bool computeelastic= true;
     5597        bool computerigid= false;
    55975598        int  glfraction=1;
    55985599
    … …  
    56405641        rho_ice=FindParam(MaterialsRhoIceEnum);
    56415642        rho_water=FindParam(MaterialsRhoSeawaterEnum);
    5642         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
     5643        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    56435644        this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
    56445645        this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
    56455646
    56465647        /*retrieve precomputed G:*/
    5647         if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
     5648        if(computerigid)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
    56485649
    56495650        /*Get area of element: precomputed in the sealevelrise_core_geometry:*/
    … …  
    56995700        _assert_(!xIsNan<IssmDouble>(bslrice));
    57005701
    5701         if(computeelastic){
     5702        if(computerigid){
    57025703                /*convert from m to kg/m^2:*/
    57035704                I=I*rho_ice*phi;
    … …  
    57255726        bool notfullygrounded=false;
    57265727        bool scaleoceanarea= false;
    5727         bool computeelastic= true;
     5728        bool computeelastic= false;
    57285729
    57295730        /*elastic green function:*/
    … …  
    57995800        IssmDouble BP;  //change in bottom pressure (Farrel and Clarke, Equ. 4)
    58005801        IssmDouble constant;
    5801         bool computeelastic= true;
     5802        bool computeelastic= false;
    58025803
    58035804        /*elastic green function:*/
    … …  
    59015902
    59025903        /*computational flags:*/
    5903         bool computeelastic= true;
     5904        bool computeelastic= false;
    59045905
    59055906        /*early return if we are not on the ocean or on an ice cap:*/

  • issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp

    r25761 r25767  
    1212#include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h"
    1313#include "../InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h"
    14                        
     14
    1515void  InputUpdateFromDakotax(FemModel* femmodel,double* variables,char* *variables_descriptors,int numdakotavariables){ /*{{{*/
    1616
    … …  
    3434
    3535        /*retrieve parameters: */
    36         femmodel->parameters->FindParam(&variable_partitions,&variable_partitions_num,NULL,NULL,QmuVariablePartitionsEnum); 
    37         femmodel->parameters->FindParam(&variable_partitions_npart,NULL,NULL,QmuVariablePartitionsNpartEnum); 
    38         femmodel->parameters->FindParam(&variable_partitions_nt,NULL,NULL,QmuVariablePartitionsNtEnum); 
    39        
    40 
    41         /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and 
    42          * for each descriptor, take the variable value and plug it into the inputs (more or less :)): 
    43          * We also start with distributed and standard values , as they tend to be used to pluck data from a multi-modle ensemble (mme) 
    44          * which can then be scaled. Doing the scaling first would be impractical, as the entire mme would have to be scaled, 
     36        femmodel->parameters->FindParam(&variable_partitions,&variable_partitions_num,NULL,NULL,QmuVariablePartitionsEnum);
     37        femmodel->parameters->FindParam(&variable_partitions_npart,NULL,NULL,QmuVariablePartitionsNpartEnum);
     38        femmodel->parameters->FindParam(&variable_partitions_nt,NULL,NULL,QmuVariablePartitionsNtEnum);
     39
     40
     41        /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and
     42         * for each descriptor, take the variable value and plug it into the inputs (more or less :)):
     43         * We also start with distributed and standard values , as they tend to be used to pluck data from a multi-modle ensemble (mme)
     44         * which can then be scaled. Doing the scaling first would be impractical, as the entire mme would have to be scaled,
    4545         * which is a waste of time:*/
    4646
    4747        variablecount=0;
    48         for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions. 
     48        for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions.
    4949
    5050                descriptor=variables_descriptors[i];
    51        
     51
    5252
    5353                /*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */
    54                 if (strncmp(descriptor,"scaled_",7)==0){ 
     54                if (strncmp(descriptor,"scaled_",7)==0){
    5555                        /*we are skipping these for now.*/
    5656                        npart=variable_partitions_npart[variablecount];
    5757                        nt=variable_partitions_nt[variablecount];
    58                                
     58
    5959                        /*increment i to skip the distributed values just collected: */
    6060                        i+=npart*nt-1; //careful, the for loop will add 1.
    … …  
    6666                        /*we are skipping these for now.*/
    6767                }
    68                
     68
    6969                else if (strncmp(descriptor,"distributed_",12)==0){
    7070                        if (VerboseQmu())_printf0_("   updating variable " << descriptor << "\n");
    71                        
     71
    7272                        /*recover partition vector: */
    7373                        variable_partition=variable_partitions[variablecount];
    7474                        npart=variable_partitions_npart[variablecount];
    7575
    76                         /*Variable is distributed. Determine root name of variable (ex: distributed_DragCoefficient_1 -> DragCoefficient). 
     76                        /*Variable is distributed. Determine root name of variable (ex: distributed_DragCoefficient_1 -> DragCoefficient).
    7777                         * Allocate distributed_values and fill the distributed_values with the next npart variables: */
    7878
    … …  
    8686
    8787                        //for (int j=0;j<npart;j++)_printf_(j << ":" << distributed_values[j] << "\n");
    88                        
     88
    8989                        //Call specialty code:
    9090                        InputUpdateSpecialtyCode(femmodel,distributed_values,variable_partition,npart,root);
    91                        
     91
    9292                        /*increment i to skip the distributed values just collected: */
    9393                        i+=npart-1; //careful, the for loop will add 1.
    94                        
     94
    9595                        /*Free allocations: */
    9696                        xDelete<double>(parameter);
    … …  
    104104                variablecount++;
    105105        }
    106        
     106
    107107        variablecount=0;
    108108        /*now deal with scaled variabes:*/
    109         for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions. 
     109        for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions.
    110110
    111111                descriptor=variables_descriptors[i];
    112        
     112
    113113                /*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */
    114114                if (strncmp(descriptor,"scaled_",7)==0){
    115                
     115
    116116                        if (VerboseQmu())_printf0_("   updating variable " << descriptor << "\n");
    117                
     117
    118118                        /*recover partition vector: */
    119119                        variable_partition=variable_partitions[variablecount];
    … …  
    121121                        nt=variable_partitions_nt[variablecount];
    122122
    123                         /* Variable is scaled, determine its root name (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the 
     123                        /* Variable is scaled, determine its root name (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the
    124124                         * distributed_values with the next npart variables coming from Dakota: */
    125125                        memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char));
    … …  
    136136                        /*increment i to skip the distributed values just collected: */
    137137                        i+=npart*nt-1; //careful, the for loop will add 1.
    138                        
     138
    139139                        /*Free allocations: */
    140140                        xDelete<double>(parameter);
    … …  
    202202        transientinput2 = femmodel->inputs->GetTransientInput(DummyEnum);
    203203
    204         for (int p=npart;p>=0;p--){
    205                 int pp=p;
    206                 if (p==npart)pp=-1; /*so, the logic is, we want to do the -1 partition first, then
    207                                                          go from npart-1 to 0 in reverse order:*/
    208                
    209                 for (int i=0;i<femmodel->elements->Size();i++){
    210                         int element_partition;
    211 
    212                         Tria*   element=xDynamicCast<Tria*>(femmodel->elements->GetObjectByOffset(i));
    213                        
    214                         element_partition= (int)variable_partition[element->Sid()];
    215                         if(element_partition!=pp)continue;
    216 
    217                         if(element_partition==-1)id=0; //grab background field
    218                         else id=distributed_values[element_partition]-1; //grab partition field
    219 
    220                         /*recover the right field from the mme: */
    221                         transientinput = datasetinput->GetTransientInputByOffset(id); _assert_(transientinput);
    222 
    223                         /*copy values from the transientinput to the final transientinput2: */
    224                         for (int j=0;j<N;j++){
    225                                 TriaInput* tria_input=transientinput->GetTriaInput(j);
    226                                 element->InputServe(tria_input);
    227                                 if(interpolationenum==P0Enum){
    228                                         value=tria_input->element_values[0];
    229                                         transientinput2->AddTriaTimeInput( j,1,&(element->lid),&value,P0Enum);
    230                                 }
    231                                 else if(interpolationenum==P1Enum){
    232 
    233                                         /*Get values and lid list*/
    234                                         const int   numvertices     = element->GetNumberOfVertices();
    235                                         int        *vertexlids      = xNew<int>(numvertices);
    236                                         int        *vertexsids      = xNew<int>(numvertices);
    237 
    238                                         /*Recover vertices ids needed to initialize inputs*/
    239                                         element->GetVerticesLidList(&vertexlids[0]);
    240                                         element->GetVerticesSidList(&vertexsids[0]);
    241                                         values=tria_input->element_values;
    242                                         transientinput2->AddTriaTimeInput( j,numvertices,vertexlids,values,P1Enum);
    243                                 }
     204        for(Object* & object : femmodel->elements->objects){
     205                Tria*   element=xDynamicCast<Tria*>(object);
     206
     207                if((int)variable_partition[element->Sid()]==-1)id=0; //grab background field
     208                else id=distributed_values[(int)variable_partition[element->Sid()]]-1; //grab partition field
     209
     210                /*recover the right field from the mme: */
     211                transientinput = datasetinput->GetTransientInputByOffset(id); _assert_(transientinput);
     212
     213                /*copy values from the transientinput to the final transientinput2: */
     214                for (int j=0;j<N;j++){
     215                        TriaInput* tria_input=transientinput->GetTriaInput(j);
     216                        element->InputServe(tria_input);
     217                        if(interpolationenum==P0Enum){
     218                                value=tria_input->element_values[0];
     219                                transientinput2->AddTriaTimeInput( j,1,&(element->lid),&value,P0Enum);
     220                        }
     221                        else if(interpolationenum==P1Enum){
     222
     223                                /*Get values and lid list*/
     224                                const int   numvertices     = element->GetNumberOfVertices();
     225                                int        *vertexlids      = xNew<int>(numvertices);
     226                                int        *vertexsids      = xNew<int>(numvertices);
     227
     228                                /*Recover vertices ids needed to initialize inputs*/
     229                                element->GetVerticesLidList(&vertexlids[0]);
     230                                element->GetVerticesSidList(&vertexsids[0]);
     231                                values=tria_input->element_values;
     232                                transientinput2->AddTriaTimeInput( j,numvertices,vertexlids,values,P1Enum);
    244233                        }
    245234                }
    … …  
    258247
    259248        /*Go through elements, copy input name to dummy, and scale it using the distributed_values and the partition vector:*/
    260         for(int i=0;i<femmodel->elements->Size();i++){
    261                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     249        for(Object* & object : femmodel->elements->objects){
     250                Element* element=xDynamicCast<Element*>(object);
    262251                element->InputScaleFromDakota(distributed_values,partition,npart,nt,name);
    263252        }
    264253
    265         /*We created a dummy input, which was a scaled copy of the name input. Now wipe 
     254        /*We created a dummy input, which was a scaled copy of the name input. Now wipe
    266255         * out the name input with the new input:*/
    267256        femmodel->inputs->ChangeEnum(DummyEnum,name);

  • issm/trunk-jpl/src/m/classes/model.m

    r25758 r25767  
    195195                        switch nargin
    196196                                case 0
    197                                         md=setdefaultparameters(md);
    198                                 case 1
    199                                         error('model constructor not supported yet');
    200 
     197                                        md=setdefaultparameters(md,'earth');
    201198                                otherwise
    202                                         error('model constructor error message: 0 of 1 argument only in input.');
    203                                 end
    204 
     199                                        options=pairoptions(varargin{:});
     200                                        planet=getfieldvalue(options,'planet','earth');
     201                                        md=setdefaultparameters(md,planet);
     202                                end
     203
     204                end
     205                %}}}
     206                function md = setdefaultparameters(md,planet) % {{{
     207
     208                        %initialize subclasses
     209                        md.mesh             = mesh2d();
     210                        md.mask             = mask();
     211                        md.constants        = constants();
     212                        md.geometry         = geometry();
     213                        md.initialization   = initialization();
     214                        md.smb              = SMBforcing();
     215                        md.basalforcings    = basalforcings();
     216                        md.friction         = friction();
     217                        md.rifts            = rifts();
     218                        md.solidearth       = solidearth(planet);
     219                        md.dsl              = dsl();
     220                        md.timestepping     = timestepping();
     221                        md.groundingline    = groundingline();
     222                        md.materials        = matice();
     223                        md.damage           = damage();
     224                        md.flowequation     = flowequation();
     225                        md.debug            = debug();
     226                        md.verbose          = verbose();
     227                        md.settings         = issmsettings();
     228                        md.toolkits         = toolkits();
     229                        md.cluster          = generic();
     230                        md.balancethickness = balancethickness();
     231                        md.stressbalance    = stressbalance();
     232                        md.hydrology        = hydrologyshreve();
     233                        md.masstransport    = masstransport();
     234                        md.thermal          = thermal();
     235                        md.steadystate      = steadystate();
     236                        md.transient        = transient();
     237                        md.levelset         = levelset();
     238                        md.calving          = calving();
     239                        md.frontalforcings  = frontalforcings();
     240                        md.gia              = giamme();
     241                        md.love             = fourierlove();
     242                        md.esa              = esa();
     243                        md.autodiff         = autodiff();
     244                        md.inversion        = inversion();
     245                        md.qmu              = qmu();
     246                        md.amr              = amr();
     247                        md.radaroverlay     = radaroverlay();
     248                        md.results          = struct();
     249                        md.outputdefinition = outputdefinition();
     250                        md.miscellaneous    = miscellaneous();
     251                        md.private          = private();
    205252                end
    206253                %}}}
    … …  
    12881335                        end
    12891336                end% }}}
    1290                 function md = setdefaultparameters(md) % {{{
    1291 
    1292                         %initialize subclasses
    1293                         md.mesh             = mesh2d();
    1294                         md.mask             = mask();
    1295                         md.constants        = constants();
    1296                         md.geometry         = geometry();
    1297                         md.initialization   = initialization();
    1298                         md.smb              = SMBforcing();
    1299                         md.basalforcings    = basalforcings();
    1300                         md.friction         = friction();
    1301                         md.rifts            = rifts();
    1302                         md.solidearth       = solidearth();
    1303                         md.dsl              = dsl();
    1304                         md.timestepping     = timestepping();
    1305                         md.groundingline    = groundingline();
    1306                         md.materials        = matice();
    1307                         md.damage           = damage();
    1308                         md.flowequation     = flowequation();
    1309                         md.debug            = debug();
    1310                         md.verbose          = verbose();
    1311                         md.settings         = issmsettings();
    1312                         md.toolkits         = toolkits();
    1313                         md.cluster          = generic();
    1314                         md.balancethickness = balancethickness();
    1315                         md.stressbalance    = stressbalance();
    1316                         md.hydrology        = hydrologyshreve();
    1317                         md.masstransport    = masstransport();
    1318                         md.thermal          = thermal();
    1319                         md.steadystate      = steadystate();
    1320                         md.transient        = transient();
    1321                         md.levelset         = levelset();
    1322                         md.calving          = calving();
    1323                         md.frontalforcings  = frontalforcings();
    1324                         md.gia              = giamme();
    1325                         md.love             = fourierlove();
    1326                         md.esa              = esa();
    1327                         md.autodiff         = autodiff();
    1328                         md.inversion        = inversion();
    1329                         md.qmu              = qmu();
    1330                         md.amr              = amr();
    1331                         md.radaroverlay     = radaroverlay();
    1332                         md.results          = struct();
    1333                         md.outputdefinition = outputdefinition();
    1334                         md.miscellaneous    = miscellaneous();
    1335                         md.private          = private();
    1336                 end
    1337                 %}}}
    13381337                function md = tetras(md,varargin) % {{{
    13391338                        %TETRAS - split 3d prismatic mesh into 3 tetrahedrons

  • issm/trunk-jpl/src/m/classes/model.py

    r25758 r25767  
    33import copy
    44import sys
     5from pairoptions import pairoptions
    56from mesh2d import mesh2d
    67from mesh3dprisms import mesh3dprisms
    … …  
    7980class model(object):
    8081    #properties
    81     def __init__(self): #{{{
    82         self.mesh = mesh2d()
    83         self.mask = mask()
    84         self.constants = constants()
    85         self.geometry = geometry()
    86         self.initialization = initialization()
    87         self.smb = SMBforcing()
    88         self.basalforcings = basalforcings()
    89         self.friction = friction()
    90         self.rifts = rifts()
    91         self.solidearth = solidearth()
    92         self.dsl = dsl()
    93         self.timestepping = timestepping()
    94         self.groundingline = groundingline()
    95         self.materials = matice()
    96         self.damage = damage()
    97         self.flowequation = flowequation()
    98         self.debug = debug()
    99         self.verbose = verbose()
    100         self.settings = issmsettings()
    101         self.toolkits = toolkits()
    102         self.cluster = generic()
    103         self.balancethickness = balancethickness()
    104         self.stressbalance = stressbalance()
    105         self.hydrology = hydrologyshreve()
    106         self.masstransport = masstransport()
    107         self.thermal = thermal()
    108         self.steadystate = steadystate()
    109         self.transient = transient()
    110         self.levelset = levelset()
    111         self.calving = calving()
    112         self.frontalforcings = frontalforcings()
    113         self.gia = giamme()
    114         self.love = fourierlove()
    115         self.esa = esa()
    116         self.autodiff = autodiff()
    117         self.inversion = inversion()
    118         self.qmu = qmu()
    119         self.amr = amr()
    120         self.radaroverlay = radaroverlay()
    121         self.results = results()
    122         self.outputdefinition = outputdefinition()
    123         self.miscellaneous = miscellaneous()
    124         self.private = private()
    125     #}}}
     82    def __init__(self, *args): #{{{
     83        self.mesh = None
     84        self.mask = None
     85        self.constants = None
     86        self.geometry = None
     87        self.initialization = None
     88        self.smb = None
     89        self.basalforcings = None
     90        self.friction = None
     91        self.rifts = None
     92        self.solidearth = None
     93        self.dsl = None
     94        self.timestepping = None
     95        self.groundingline = None
     96        self.materials = None
     97        self.damage = None
     98        self.flowequation = None
     99        self.debug = None
     100        self.verbose = None
     101        self.settings = None
     102        self.toolkits = None
     103        self.cluster = None
     104        self.balancethickness = None
     105        self.stressbalance = None
     106        self.hydrology = None
     107        self.masstransport = None
     108        self.thermal = None
     109        self.steadystate = None
     110        self.transient = None
     111        self.levelset = None
     112        self.calving = None
     113        self.frontalforcings = None
     114        self.gia = None
     115        self.love = None
     116        self.esa = None
     117        self.autodiff = None
     118        self.inversion = None
     119        self.qmu = None
     120        self.amr = None
     121        self.radaroverlay = None
     122        self.results = None
     123        self.outputdefinition = None
     124        self.miscellaneous = None
     125        self.private = None
     126
     127        nargs = len(args)
     128
     129        if nargs == 0:
     130            self.setdefaultparameters('earth')
     131        else:
     132            self.setdefaultparameters(args[0])
     133            options = pairoptions(*args)
     134            planet = options.getfieldvalue('planet', 'earth')
     135            self.setdefaultparameters(planet)
     136#}}}
    126137
    127138    def properties(self):  # {{{
    … …  
    224235    # }}}
    225236
     237    def setdefaultparameters(self, planet): #{{{
     238        self.mesh = mesh2d()
     239        self.mask = mask()
     240        self.constants = constants()
     241        self.geometry = geometry()
     242        self.initialization = initialization()
     243        self.smb = SMBforcing()
     244        self.basalforcings = basalforcings()
     245        self.friction = friction()
     246        self.rifts = rifts()
     247        self.solidearth = solidearth(planet)
     248        self.dsl = dsl()
     249        self.timestepping = timestepping()
     250        self.groundingline = groundingline()
     251        self.materials = matice()
     252        self.damage = damage()
     253        self.flowequation = flowequation()
     254        self.debug = debug()
     255        self.verbose = verbose()
     256        self.settings = issmsettings()
     257        self.toolkits = toolkits()
     258        self.cluster = generic()
     259        self.balancethickness = balancethickness()
     260        self.stressbalance = stressbalance()
     261        self.hydrology = hydrologyshreve()
     262        self.masstransport = masstransport()
     263        self.thermal = thermal()
     264        self.steadystate = steadystate()
     265        self.transient = transient()
     266        self.levelset = levelset()
     267        self.calving = calving()
     268        self.frontalforcings = frontalforcings()
     269        self.gia = giamme()
     270        self.love = fourierlove()
     271        self.esa = esa()
     272        self.autodiff = autodiff()
     273        self.inversion = inversion()
     274        self.qmu = qmu()
     275        self.amr = amr()
     276        self.radaroverlay = radaroverlay()
     277        self.results = results()
     278        self.outputdefinition = outputdefinition()
     279        self.miscellaneous = miscellaneous()
     280        self.private = private()
     281    #}}}
     282
    226283    def checkmessage(self, string):  # {{{
    227284        print("model not consistent: {}".format(string))

  • issm/trunk-jpl/src/m/classes/solidearth.m

    r25763 r25767  
    66classdef solidearth
    77        properties (SetAccess=public)
    8                 initialsealevel               = NaN;
     8                initialsealevel        = NaN;
    99                settings               = solidearthsettings();
    1010                external               = [];
    … …  
    2222                        switch nargin
    2323                                case 0
    24                                         self=setdefaultparameters(self);
     24                                        self=setdefaultparameters(self,'earth');
     25                                case 1
     26                                        self=setdefaultparameters(self,varargin{:});
    2527                                otherwise
    26                                         error('constructor not supported');
     28                                        error('solidearth constructor error message: zero or one argument only!');
    2729                        end
    2830                end % }}}
    29                 function self = setdefaultparameters(self) % {{{
    30                
    31                 %output default:
    32                 self.requested_outputs={'default'};
     31                function self = setdefaultparameters(self,planet) % {{{
    3332
    34                 %transitions should be a cell array of vectors:
    35                 self.transitions={};
    36                
    37                 %no partitions requested for barystatic contribution:
    38                 self.partitionice=[];
    39                 self.partitionhydro=[];
     33                        %output default:
     34                        self.requested_outputs={'default'};
    4035
    41                 %no external solutions by defalt:
    42                 self.external=[];
     36                        %transitions should be a cell array of vectors:
     37                        self.transitions={};
     38                       
     39                        %no partitions requested for barystatic contribution:
     40                        self.partitionice=[];
     41                        self.partitionhydro=[];
    4342
    44                 %earth radius
    45                 self.planetradius= planetradius('earth');
     43                        %no external solutions by defalt:
     44                        self.external=[];
     45
     46                        %earth radius
     47                        self.planetradius= planetradius(planet);
    4648               
    4749                end % }}}

  • issm/trunk-jpl/src/m/classes/solidearth.py

    r25764 r25767  
    3232        self.partitionhydro     = []
    3333
    34         nargin = len(args)
     34        nargs = len(args)
    3535
    36         if nargin == 0:
    37             self.setdefaultparameters()
     36        if nargs == 0:
     37            self.setdefaultparameters('earth')
     38        elif nargs == 1:
     39            self.setdefaultparameters(args[0])
    3840        else:
    39             raise Exception('constructor not supported')
     41            raise Exception('solidearth constructor error message: zero or one argument only!')
    4042    #}}}
    4143
    … …  
    5961    #}}}
    6062
    61     def setdefaultparameters(self):  # {{{
     63    def setdefaultparameters(self, planet):  # {{{
    6264        # Default output
    6365        self.requested_outputs = ['default']
    … …  
    7476
    7577        # Earth radius
    76         self.planetradius = planetradius('earth')
     78        self.planetradius = planetradius(planet)
    7779    #}}}
    7880

  • issm/trunk-jpl/src/m/classes/solidearthsettings.m

    r25763 r25767  
    7171                        md = checkfield(md,'fieldname','solidearth.settings.glfraction','values',[0 1]);
    7272
     73                        %checks on computational flags
     74                        if self.elastic==1 & self.rigid==0,
     75                                error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set');
     76                        end
     77
    7378                        %a coupler to a planet model is provided.
    7479                        if self.isgrd,

  • issm/trunk-jpl/src/m/classes/solidearthsettings.py

    r25766 r25767  
    9292        md = checkfield(md, 'fieldname', 'solidearth.settings.glfraction', 'values', [0, 1])
    9393
     94        # Checks on computational flags
     95        if self.elastic and not self.rigid:
     96            raise Exception('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set')
     97
    9498        # A coupler to planet model is provided
    9599        if self.isgrd:

  • issm/trunk-jpl/src/m/geometry/planetradius.py

    r24902 r25767  
    1111
    1212    if planet == 'earth':
    13         radius = 6.371012 * pow(10, 6)
     13        radius = 6.371012e6
     14    elif planet == 'europa':
     15        radius = 1.5008e6
    1416    else:
    1517        raise TypeError("planet type %s not supported yet!" % planet)
Note: See TracChangeset for help on using the changeset viewer.