Changeset 26468 for issm/trunk-jpl/src


Ignore:
Timestamp:
10/04/21 06:35:33 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: removing extraneous lines

Location:
issm/trunk-jpl/src/c
Files:
59 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

    r26362 r26468  
    150150}/*}}}*/
    151151ElementMatrix* AdjointHorizAnalysis::CreateKMatrixHO(Element* element){/*{{{*/
    152        
     152
    153153        /* Check if ice in element */
    154154        if(!element->IsIceInElement()) return NULL;

  • TabularUnified issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

    r26223 r26468  
    7070                element->SetElementInput(inputs,DamageFEnum,0.);
    7171        }
    72 
    7372
    7473        /*What input do I need to run my damage evolution model?*/

  • TabularUnified issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

    r26090 r26468  
    535535        int iseplthickcomp;
    536536
    537 
    538537        /*Skip if we don't want to compute thicknesses*/
    539538        femmodel->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum);
    … …  
    558557        femmodel->parameters->FindParam(&gravity,ConstantsGEnum);
    559558        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    560 
    561559
    562560        for(Object* & object : femmodel->elements->objects){

  • TabularUnified issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp

    r26463 r26468  
    562562        /*Skip if water or ice shelf element*/
    563563        if(element->IsAllFloating()) return;
    564        
     564
    565565        /*Intermediary*/
    566566        IssmDouble phi_0, phi_m, p_i;

  • TabularUnified issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp

    r26329 r26468  
    307307        /*Add input to the element: */
    308308        element->AddInput(WatercolumnEnum,watercolumn,element->GetElementType());
    309        
     309
    310310        /*Also update the hydrological loads for the sealevel core: */
    311311        IssmDouble* oldwatercolumn      = xNew<IssmDouble>(numnodes);

  • TabularUnified issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp

    r26075 r26468  
    4242        iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",WatercolumnEnum);
    4343
    44 
    4544}/*}}}*/
    4645void HydrologyTwsAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    … …  
    8584        int numnodes = element->GetNumberOfNodes();
    8685        IssmDouble* watercolumn = xNew<IssmDouble>(numnodes);
    87        
     86
    8887        /*Use the dof list to index into the solution vector: */
    8988        for(int i=0;i<numnodes;i++){

  • TabularUnified issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

    r26465 r26468  
    576576                IssmDouble dt                 = femmodel->parameters->FindParam(TimesteppingTimeStepEnum);
    577577
    578 
    579578                for(Object* & object : femmodel->elements->objects){
    580579                        Element* element   = xDynamicCast<Element*>(object);

  • TabularUnified issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r26090 r26468  
    234234        if(grdmodel==IvinsEnum) inputs->SetTransientInput(TransientAccumulatedDeltaIceThicknessEnum,NULL,0);
    235235
    236 
    237236}/*}}}*/
    238237void MasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    … …  
    816815                element->FindParam(&count,SealevelchangeRunCountEnum);
    817816        }
    818                
     817
    819818        /*Fetch dof list and allocate solution vector*/
    820819        int *doflist = NULL;
    … …  
    912911        element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
    913912        element->AddBasalInput(BaseEnum,newbase,P1Enum);
    914        
    915913
    916914        /*Free ressources:*/

  • TabularUnified issm/trunk-jpl/src/c/analyses/OceantransportAnalysis.cpp

    r26099 r26468  
    6262        iomodel->FetchDataToInput(inputs,elements,"md.initialization.str",StrEnum);
    6363
    64 
    6564}/*}}}*/
    6665void OceantransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    … …  
    7574                IssmDouble modelid;
    7675                int nummodels;
    77                
     76
    7877                /*create double param, not int param, because Dakota will be updating it as
    7978                 * a double potentially: */

  • TabularUnified issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp

    r26296 r26468  
    4949                iomodel->FetchDataToInput(inputs,elements,"md.solidearth.external.geoid",SolidearthExternalGeoidRateEnum);
    5050
    51 
    5251                /*Resolve Mmes using the modelid, if necessary:*/
    5352                if (inputs->GetInputObjectEnum(SolidearthExternalDisplacementEastRateEnum)==DatasetInputEnum){
    5453                        int modelid;
    55                        
     54
    5655                        /*retrieve model id: */
    5756                        iomodel->FetchData(&modelid,"md.solidearth.external.modelid");
    58                
     57
    5958                        /*replace dataset of forcings with only one, the modelid'th:*/
    6059                        MmeToInputFromIdx(inputs,elements,modelid,SolidearthExternalDisplacementNorthRateEnum, P1Enum);
    … …  
    7473        iomodel->ConstantToInput(inputs,elements,0.,DeltaIceThicknessEnum,P1Enum);
    7574        iomodel->ConstantToInput(inputs,elements,0.,DeltaBottomPressureEnum,P1Enum);
    76 
    7775
    7876}/*}}}*/
    … …  
    182180                xDelete<IssmDouble>(partitionice);
    183181        }
    184        
     182
    185183        iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro");
    186184        parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro));
    … …  
    208206        BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel);
    209207        parameters->AddObject(new GenericParam<BarystaticContributions*>(barystaticcontributions,BarystaticContributionsEnum));
    210        
     208
    211209        /*Deal with external multi-model ensembles: {{{*/
    212210        if(isexternal){
    … …  
    425423                        xDelete<int>(displs);
    426424
    427                        
    428 
    429425                        /*Avoid singularity at 0: */
    430426                        G_gravi[0]=G_gravi[1];
    … …  
    436432                                }
    437433                        }
    438 
    439                        
    440434
    441435                        /*Reinterpolate viscoelastic green kernels onto a regular gridded time
    … …  
    484478                                        }
    485479                                }
    486                                                
     480
    487481                        }
    488482                        else if(elastic){

  • TabularUnified issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r26208 r26468  
    232232        parameters->AddObject(new IntParam(SmbStepsPerStepEnum,smbslices));
    233233
    234 
    235234        parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));
    236235

  • TabularUnified issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r26463 r26468  
    11091109                #endif
    11101110
    1111 
    11121111                if(is_schur_cg_solver)
    11131112                 solutionsequence_schurcg(femmodel);
    … …  
    20962095        /*Use the dof list to index into the solution vector: */
    20972096        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    2098        
     2097
    20992098        /*Transform solution in Cartesian Space*/
    21002099        if(dim==2) basalelement->TransformSolutionCoord(&values[0],XYEnum);
    … …  
    21342133        }
    21352134        element->AddBasalInput(VelEnum,vel,element->GetElementType());
    2136        
     2135
    21372136        /*Free ressources:*/
    21382137        xDelete<IssmDouble>(vel);
    … …  
    27542753/*MLHO*/
    27552754ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/
    2756        
     2755
    27572756        //_error_("Mono Layer Higher-Order called, not fully tested. If you are sure of using it, comment this line.");
    2758        
     2757
    27592758        /* Check if ice in element */
    27602759        if(!element->IsIceInElement()) return NULL;
    … …  
    28182817                gauss = element->NewGauss(2);
    28192818        }
    2820        
     2819
    28212820        /* Start  looping on the number of gaussian points: */
    28222821        while(gauss->next()){
    … …  
    28502849        xDelete<IssmDouble>(basis);
    28512850        return Ke;
    2852 
    28532851
    28542852        //itapopo OLD - testing above
    … …  
    28732871        //ElementMatrix* Ke    = basalelement->NewElementMatrix(MLHOApproximationEnum);
    28742872        ElementMatrix* KeSSA = CreateKMatrixSSAFriction(basalelement); //only to get K11 and K33
    2875        
     2873
    28762874        /*Fetch number of nodes and dof for this finite element*/
    28772875        //int numnodes = basalelement->GetNumberOfNodes();
    … …  
    29692967                     2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
    29702968                     )*(n+1)/(n+2);//K14
    2971                                
     2969
    29722970                                Ke->values[(4*i+1)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity[1]*thickness*(
    29732971                                                        4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
    … …  
    29852983                                                        2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
    29862984                                                        )*2*pow(n+1,2)/((2*n+3)*(n+2));//K24
    2987                                
     2985
    29882986                                Ke->values[(4*i+2)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity[0]*thickness*(
    29892987                                                        2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
    … …  
    30173015                }
    30183016        }
    3019  
     3017
    30203018        /*Transform Coordinate System*/
    30213019        //basalelement->TransformStiffnessMatrixCoord(Ke,XYMLHOEnum);
    … …  
    30513049                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    30523050        }
    3053        
     3051
    30543052        /*compute all load vectors for this element*/
    30553053        ElementVector* pe1=CreatePVectorMLHODrivingStress(basalelement);
    … …  
    30923090                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
    30933091                n_input->GetInputValue(&n,gauss);
    3094                
     3092
    30953093                for(int i=0;i<numnodes;i++){// per node: vx (basal vx), vshx, vy (basal vy), vshy
    30963094                        pe->values[i*4+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F1
    … …  
    31503148                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
    31513149                element->NodalFunctions(basis,gauss);
    3152                
     3150
    31533151                b = base-sealevel; // ice base shifted by the sea level
    31543152                s = thickness+b;   // ice surface shifted by the sea level
    … …  
    32493247        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
    32503248        IssmDouble* n                    = xNew<IssmDouble>(numnodes);
    3251        
     3249
    32523250        /*Use the dof list to index into the solution vector: */
    32533251        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    … …  
    32833281        element->AddBasalInput(VxShearEnum,vshx,element->GetElementType());
    32843282        element->AddBasalInput(VyShearEnum,vshy,element->GetElementType());
    3285        
     3283
    32863284        /*Add vx and vy as inputs to the tria element (base velocities): */
    32873285        element->AddBasalInput(VxBaseEnum,vbx,element->GetElementType());
    32883286        element->AddBasalInput(VyBaseEnum,vby,element->GetElementType());
    3289        
     3287
    32903288        /*Add vx and vy as inputs to the tria element (surface velocities): */
    32913289        element->AddBasalInput(VxSurfaceEnum,vsx,element->GetElementType());
    … …  
    32983296                vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2);
    32993297        }
    3300                
     3298
    33013299        /*Get Vz and compute vel (vertically averaged vel)*/
    33023300        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
    … …  
    33723370
    33733371        solution->SetValues(numdof,doflist,values,INS_VAL);
    3374        
     3372
    33753373        /*Free ressources:*/
    33763374        delete gauss;
    … …  
    44064404        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
    44074405
    4408 
    44094406        /* Start  looping on the number of gaussian points: */
    44104407        Gauss* gauss = element->NewGauss(5);
    … …  
    44144411                element->NodalFunctionsPressure(pbasis,gauss);
    44154412                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    4416 
    44174413
    44184414                if(dim==3 || dim==2){
    … …  
    54855481                element->NodalFunctionsPressure(pbasis,gauss);
    54865482                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    5487                
     5483
    54885484                for(int i=0;i<vnumnodes;i++){
    54895485                        for(int j=0;j<vnumnodes;j++){

  • TabularUnified issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

    r26047 r26468  
    117117        iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
    118118        #endif
    119 
    120119
    121120        /*Add basal forcings to compute melt rate*/
    … …  
    583582        }
    584583
    585 
    586584        #ifdef INWOOVZ
    587585        IssmDouble*  thickness = xNew<IssmDouble>(numnodes);
    … …  
    624622        xDelete<IssmDouble>(thickness);
    625623        #endif
    626 
    627624
    628625        /*Now Compute vel*/

  • TabularUnified issm/trunk-jpl/src/c/classes/AmrBamg.cpp

    r23501 r26468  
    8686/*Methods*/
    8787void AmrBamg::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/
    88        
     88
    8989        /*Delete previous mesh and keep the entire mesh*/
    9090        if(this->elementslist) xDelete<int>(this->elementslist);
    … …  
    9999}/*}}}*/
    100100void AmrBamg::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/
    101        
     101
    102102        /*Get the entire mesh*/
    103103        *elementslist_out               = this->elementslist;
    … …  
    128128}/*}}}*/
    129129void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int** pdatalist,IssmDouble** pxylist,int** pelementslist){/*{{{*/
    130        
     130
    131131        /*Intermediaries*/
    132132        BamgGeom* geomout=new BamgGeom();
    … …  
    142142        this->options->field                     = field;
    143143        this->options->hmaxVertices = hmaxVertices;
    144        
     144
    145145        /*remesh*/
    146146        if(this->previousmesh){

  • TabularUnified issm/trunk-jpl/src/c/classes/AmrNeopz.cpp

    r25454 r26468  
    137137/*Mesh refinement methods*/
    138138void AmrNeopz::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/
    139    
     139
    140140   /*Delete previous mesh and keep the entire mesh*/
    141141   if(this->elementslist) xDelete<int>(this->elementslist);
    … …  
    150150}/*}}}*/
    151151void AmrNeopz::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/
    152    
     152
    153153   /*Get the entire mesh*/
    154154   *elementslist_out    = this->elementslist;
    … …  
    394394                if(type>1) break;
    395395        }
    396        
     396
    397397        /*Verify and return*/
    398398        if(type>1) type=2;
    … …  
    416416
    417417        count=1;
    418        
     418
    419419        while(count>0){
    420420                count=0;
    … …  
    438438                                }
    439439                        }
    440                        
     440
    441441                        /*refine the element if requested*/
    442442                        if(type==2){gmesh->Element(i)->Divide(sons);    count++;}
    … …  
    926926        this->fathermesh = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::ReadFromFile());
    927927        TPZPersistenceManager::CloseRead();
    928        
     928
    929929        TPZPersistenceManager::OpenRead(previousmeshfile);
    930930        this->previousmesh = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::ReadFromFile());
    … …  
    933933        if(!amrifstream.is_open()) _error_("amr ifstream is not open!");
    934934        amrifstream.seekg(0);
    935        
     935
    936936        this->sid2index.clear();
    937937        this->index2sid.clear();
    938938        this->specialelementsindex.clear();
    939        
     939
    940940        amrifstream>>size;
    941941        for(int i=0;i<size;i++){
    … …  
    943943                this->sid2index.push_back(value);
    944944        }
    945        
     945
    946946        amrifstream>>size;
    947947        for(int i=0;i<size;i++){
    … …  
    949949                this->index2sid.push_back(value);
    950950        }
    951        
     951
    952952        amrifstream>>size;
    953953        for(int i=0;i<size;i++){
    … …  
    982982                }
    983983        }
    984        
     984
    985985        if(this->index2sid.size()>0){
    986986                amrofstream << this->index2sid.size() << std::endl;

  • TabularUnified issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp

    r26287 r26468  
    6666
    6767        IssmDouble  sumice,sumhydro,sumocean;
    68        
     68
    6969        ice->Assemble();
    7070        hydro->Assemble();
    … …  
    9090        cumocean->Sum(&sumocean);
    9191
    92 
    9392        return sumice+sumhydro+sumocean;
    9493
    … …  
    10099        cumhydro->AXPY(hydro,1);
    101100
    102 
    103101} /*}}}*/
    104102void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/
    105        
     103
    106104        int id;
    107105        if(nice){
    … …  
    128126                ocean->SetValue(0,oceanvalue,ADD_VAL);
    129127        }
    130                
    131128
    132129} /*}}}*/
    133130void BarystaticContributions::Reset(){ /*{{{*/
    134131
    135 
    136132        ice->Set(0.);
    137133        hydro->Set(0.);
    138134        ocean->Set(0.);
    139        
     135
    140136} /*}}}*/
    141137void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/

  • TabularUnified issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp

    r26090 r26468  
    9696IssmDouble Cfsurfacelogvel::Response(FemModel* femmodel){/*{{{*/
    9797
    98 
    9998        /*recover time parameters: */
    10099        IssmDouble time;

  • TabularUnified issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp

    r26222 r26468  
    124124                if(VerboseQmu()) _printf0_("compute dakota responses:\n");
    125125                femmodel->DakotaResponsesx(responses,responses_descriptors,numresponsedescriptors,numFns);
    126                
     126
    127127                /*Output results for this iteration: */
    128128                if(VerboseQmu()) _printf0_("output results for this iteration: \n");

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

    r26329 r26468  
    12861286/*}}}*/
    12871287void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type){/*{{{*/
    1288 
    12891288
    12901289        switch(type){
    … …  
    16801679                        if(M==1)times[0]=0;
    16811680                        if(M==2)for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1682                        
     1681
    16831682                        inputs->SetTransientInput(vector_enum,times,N);
    16841683                        TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
    1685                
     1684
    16861685                        for(int t=0;t<N;t++){
    16871686                                value=vector[t]; //values are on the first line, times are on the second line
    … …  
    32623261        }
    32633262
    3264 
    32653263        /*Assign output pointer*/
    32663264
    … …  
    33233321                                case P1Enum:{
    33243322                                        const int NUM_VERTICES = this->GetNumberOfVertices();
    3325 
    3326 
    33273323
    33283324                                        this->GetVerticesSidList(&sidlist[0]);

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

    r26362 r26468  
    19401940        }
    19411941
    1942 
    19431942        /*Some checks in debugging mode*/
    19441943        _assert_(s1>=0 && s1<=1.);
    … …  
    20552054                _error_("case not possible");
    20562055        }
    2057 
    20582056
    20592057        /*Some checks in debugging mode*/
    … …  
    23772375        if(start==+1 && !IsOnSurface()) return;
    23782376
    2379 
    23802377        /*Get original input*/
    23812378        Input* input = this->GetInput(enum_type);
    … …  
    27382735        IssmDouble  movingfrontvy[NUMVERTICES];
    27392736        IssmDouble  vel;
    2740        
     2737
    27412738        /* Get node coordinates and dof list: */
    27422739        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    40334030        }
    40344031
    4035 
    40364032        /*Some checks in debugging mode*/
    40374033        _assert_(s1>=0 && s1<=1.);
    … …  
    41474143                _error_("case not possible");
    41484144        }
    4149 
    41504145
    41514146        /*Some checks in debugging mode*/
    … …  
    48994894        else if(type==TransientInputEnum){
    49004895
    4901 
    49024896                IssmDouble* steps=NULL;
    49034897                int nsteps;

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

    r26464 r26468  
    323323        IssmDouble  epse_2,groundedice,bed,sealevel;            // added sealevel
    324324
    325 
    326325        /* Get node coordinates and dof list: */
    327326        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    337336        Input* n_input  = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
    338337        Input* sl_input  = this->GetInput(SealevelEnum); _assert_(sl_input);
    339 
    340 
    341338
    342339        /* Start looping on the number of vertices: */
    … …  
    413410        int         crevasse_opening_stress;
    414411
    415 
    416412        /*reset if no ice in element*/
    417413        if(!this->IsIceInElement()){
    … …  
    426422                return;
    427423        }
    428 
    429424
    430425        /*retrieve the type of crevasse_opening_stress*/
    … …  
    13211316                area=this->GetAreaSpherical();
    13221317                vareae->SetValue(this->sid,area,INS_VAL);
    1323                
     1318
    13241319                /*in addition, put in in the inputs:*/
    13251320                this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
    … …  
    13501345        vlate->SetValue(this->sid,late,INS_VAL);
    13511346        vareae->SetValue(this->sid,area,INS_VAL);
    1352                
     1347
    13531348        return;
    13541349}
    … …  
    17661761/*}}}*/
    17671762void       Tria::GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){/*{{{*/
    1768        
     1763
    17691764        /*Computeportion of the element that is grounded*/
    17701765        bool               trapezeisnegative=true; //default value
    … …  
    17771772        IssmDouble loadweights_g[NUMVERTICES];
    17781773        IssmDouble total_weight=0;
    1779        
     1774
    17801775        _assert_(!xIsNan<IssmDouble>(gl[0]));
    17811776        _assert_(!xIsNan<IssmDouble>(gl[1]));
    … …  
    18201815        if(trapezeisnegative) phi=1-f1*f2;
    18211816        else phi=f1*f2;
    1822        
    18231817
    18241818        /*Compute weights:*/
    18251819        gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
    1826        
     1820
    18271821        total_weight=0;
    18281822        for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
    … …  
    18481842/*}}}*/
    18491843IssmDouble Tria::GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){/*{{{*/
    1850        
     1844
    18511845        IssmDouble x1,y1,z1,x2,y2,z2,x3,y3,z3;
    18521846        IssmDouble arc12,arc23,arc31,semi_peri,excess;
    … …  
    19291923        IssmDouble levelset[NUMVERTICES];
    19301924        IssmDouble planetradius;
    1931        
     1925
    19321926        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
    19331927
    … …  
    19401934        Element::GetInputListOnVertices(&levelset[0],levelsetenum);
    19411935        if(flip1)for(int i=0;i<NUMVERTICES;i++)levelset[i]=-levelset[i];
    1942        
     1936
    19431937        //compute sea level load weights
    19441938        this->GetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset);
    … …  
    19511945} /*}}}*/
    19521946void       Tria::GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelset1enum, int levelset2enum){ /*{{{*/
    1953 
    19541947
    19551948        bool istrapeze1, istrapeze2;
    … …  
    20192012        }
    20202013
    2021 
    20222014        //If everyone is negative, no need to calculate any fraction
    20232015        if (levelset1[0]<=0 && levelset1[1]<=0 && levelset1[2]<=0 && levelset2[0]<=0 && levelset2[1]<=0 && levelset2[2]<=0) {
    … …  
    20282020                return;
    20292021        }
    2030        
     2022
    20312023        /*recover planet radius:*/
    20322024        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
    2033 
    20342025
    20352026        //If just one levelset is all negative, just take the partitioning of the other, no interaction between them
    … …  
    20472038        }
    20482039
    2049 
    20502040        this->GetFractionGeometry(&weights1[0],&phi1,&point1,&d,&e,&istrapeze1,levelset1);
    20512041        this->GetFractionGeometry(&weights2[0],&phi2,&point2,&f,&g,&istrapeze2,levelset2);
    2052 
    20532042
    20542043        //Early return if levelsets are not independent
    … …  
    20692058        }
    20702059
    2071        
    20722060        ::GetVerticesCoordinates(&xyz0[0][0],vertices,NUMVERTICES); // initial triangle
    20732061
    … …  
    21172105                h2= (d*(f+e-1))/(g*(e-d) +f*d);
    21182106        }
    2119 
    21202107
    21212108        //interpolant weights of each point. Any field F[0,1,2] provided at the original vertices [0,1,2] will be equal on point k to sum_i (F[i] * w[k][i])
    … …  
    22272214                                } /*}}}*/
    22282215                        }
    2229                        
     2216
    22302217                }/*}}}*/
    22312218                else{ /*{{{*/
    … …  
    24982485}/*}}}*/
    24992486Input*    Tria::GetInput(int inputenum){/*{{{*/
    2500                
     2487
    25012488        /*Get Input from dataset*/
    25022489        if(this->iscollapsed){
    … …  
    40814068        IssmDouble  movingfrontvy[NUMVERTICES];
    40824069        IssmDouble  vel;
    4083        
     4070
    40844071        /* Get node coordinates and dof list: */
    40854072        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    43324319                }
    43334320                for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
    4334                                
     4321
    43354322                movingfrontvx[iv] = w[0];
    43364323                movingfrontvy[iv] = w[1];               
    … …  
    61476134        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    61486135        this->parameters->FindParam(&yts,ConstantsYtsEnum);
    6149        
     6136
    61506137        /*recover gia solution parameters: */
    61516138        int cross_section_shape;
    … …  
    61586145        /*recover material parameters: */
    61596146        IssmDouble rho_ice                   = FindParam(MaterialsRhoIceEnum);
    6160        
     6147
    61616148        /*recover mantle and lithosphere material properties:*/
    61626149        int numlayers=litho->numlayers;
    … …  
    61796166        int         numtimes;
    61806167        this->GetInputAveragesUpToCurrentTime(TransientAccumulatedDeltaIceThicknessEnum,&hes,&times,&numtimes,currenttime);
    6181 
    61826168
    61836169        /*pull area of this Tria: */
    … …  
    62746260        IssmDouble* H_viscoelastic_precomputed=NULL;
    62756261        #endif
    6276        
     6262
    62776263        /*viscoelastic green function:*/
    62786264        int index;
    … …  
    63426328                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
    63436329                parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
    6344                
     6330
    63456331                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
    63466332                parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
    … …  
    63736359        }
    63746360
    6375 
    63766361        constant=3/rho_earth/planetarea;
    63776362
    … …  
    63936378                        IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
    63946379                        IssmDouble longe= atan2(yye[e],xxe[e]);
    6395                        
     6380
    63966381                        lati=latitude[i];
    63976382                        longi=longitude[i];
    … …  
    64046389                        _assert_(index>=0 && index<M);
    64056390
    6406 
    64076391                        lincoef=doubleindex-index; //where between index and index+1 is alpha
    64086392                        if (index==M-1){ //avoids out of bound case
    … …  
    64106394                                lincoef=1;
    64116395                        }
    6412                        
    64136396
    64146397                        if(horiz){
    … …  
    65916574        /* Classic buildup of load weights, centroids and areas *for elements which are fully inside a mask.
    65926575         * At the same time, we'll tag the elements that are fractionally only inside a mask*/
    6593        
     6576
    65946577        IssmDouble loadweights[3]={0};
    65956578        IssmDouble area;
    … …  
    66106593        /*constants:*/
    66116594        IssmDouble constant=0;
    6612        
     6595
    66136596        /*recover parameters:*/
    66146597        this->parameters->FindParam(&computeice,TransientIsmasstransportEnum);
    … …  
    66716654                slgeom->nsubel[SLGEOM_OCEAN]++;
    66726655        }
    6673                
     6656
    66746657        /*early return if we are not on an ice sheet , and we are not requesting
    66756658         *hydrology or bottom pressure loads :*/
    … …  
    67486731                }
    67496732        }
    6750        
     6733
    67516734        /*Deal with ice loads if we are on grounded ice:*/
    67526735        if(isice && !isoceanonly && computeice){
    … …  
    67916774                }
    67926775        }
    6793        
     6776
    67946777}
    67956778/*}}}*/
    … …  
    67986781        /* Classic buildup of load weights, centroids and areas *for elements which are fully inside a mask.
    67996782         * At the same time, we'll tag the elements that are fractionally only inside a mask*/
    6800        
     6783
    68016784        IssmDouble loadweights[3]={0};
    68026785        IssmDouble area,loadarea;
    … …  
    68086791        IssmDouble constant;
    68096792        IssmDouble nanconstant=NAN;
    6810        
     6793
    68116794        /*get vertex and area information:*/
    68126795        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    68216804        this->AddInput(SealevelBarystaticOceanLongbarEnum,&longbar,P0Enum);
    68226805        #endif
    6823        
     6806
    68246807        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    68256808                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    6826                
     6809
    68276810                this->GetNodalWeightsAndAreaAndCentroidsFromLeveset(&loadweightsocean[0],&loadareaocean,&latbar, &longbar, slgeom->late[this->lid], slgeom->longe[this->lid], area, MaskOceanLevelsetEnum);
    68286811                slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=loadareaocean;
    … …  
    68856868                this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum);
    68866869                this->AddInput(SealevelBarystaticHydroAreaEnum,&loadarea,P0Enum);
    6887                
     6870
    68886871                this->AddInput(SealevelBarystaticHydroLatbarEnum,&latbar,P0Enum);
    68896872                this->AddInput(SealevelBarystaticHydroLongbarEnum,&longbar,P0Enum);
    … …  
    68916874                #endif
    68926875        }
    6893        
     6876
    68946877}
    68956878/*}}}*/
    … …  
    69276910        IssmDouble** GEsubel=NULL;
    69286911        #endif
    6929        
     6912
    69306913        /*viscoelastic green function:*/
    69316914        int index;
    … …  
    70106993                }
    70116994        }
    7012        
     6995
    70136996        //Allocate:
    70146997        for(int l=0;l<SLGEOM_NUMLOADS;l++){
    … …  
    71337116/*}}}*/
    71347117void       Tria::SealevelchangeUpdateViscousFields(){ /*{{{*/
    7135        
     7118
    71367119        /*Inputs:*/
    71377120        IssmDouble* viscousRSL=NULL;
    … …  
    71487131        IssmDouble  lincoeff=0;
    71497132        int horiz;
    7150                
     7133
    71517134        this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
    7152        
     7135
    71537136        if(viscous){
    71547137                this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    7155                
     7138
    71567139                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    71577140                this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
    … …  
    71657148                        this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&dummy);
    71667149                }
    7167 
    71687150
    71697151                bool foundtime=false;
    … …  
    71947176                viscousindex=newindex+offset;
    71957177
    7196                
    71977178                this->parameters->SetParam(viscousindex,SealevelchangeViscousIndexEnum);
    71987179                this->parameters->SetParam(viscoustimes,viscousnumsteps,SealevelchangeViscousTimesEnum);
    … …  
    72017182                xDelete<IssmDouble>(viscoustimes);
    72027183        }
    7203 
    72047184
    72057185}
    … …  
    72547234        bslchydro = -slgeom->LoadArea[SLGEOM_WATER][this->lid]*Wavg;
    72557235        bslcbp =    -slgeom->LoadArea[SLGEOM_OCEAN][this->lid]*BPavg;
    7256        
     7236
    72577237        _assert_(!xIsNan<IssmDouble>(bslcice));
    72587238        _assert_(!xIsNan<IssmDouble>(bslchydro));
    72597239        _assert_(!xIsNan<IssmDouble>(bslcbp));
    7260        
     7240
    72617241        /*Plug values into subelement load vector:*/
    72627242        if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
    … …  
    72897269        IssmDouble* Gsub[SLGEOM_NUMLOADS];
    72907270        bool computefuture=false;
    7291        
     7271
    72927272        bool sal = false;
    72937273        bool viscous = false;
    … …  
    72997279        IssmDouble Grotm2[3];
    73007280        IssmDouble Grotm3[3];
    7301                
     7281
    73027282        this->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum);
    73037283        this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
    … …  
    73187298                Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
    73197299                Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
    7320                
     7300
    73217301                for(int i=0;i<NUMVERTICES;i++){
    73227302                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    … …  
    73337313        IssmDouble oceanarea=0;
    73347314        IssmDouble rho_water;
    7335        
     7315
    73367316        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    73377317
    … …  
    73527332        }
    73537333        else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    7354 
    73557334
    73567335        /*add ocean area into a global oceanareas vector:*/
    … …  
    73627341                }
    73637342        }
    7364        
     7343
    73657344        return;
    73667345} /*}}}*/
    … …  
    74097388        this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    74107389        this->parameters->FindParam(&planethasocean,SolidearthSettingsGrdOceanEnum);
    7411        
    7412                
     7390
    74137391        if(sal){
    74147392
    … …  
    74507428                Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
    74517429                Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
    7452                
     7430
    74537431                for(int i=0;i<NUMVERTICES;i++){
    74547432                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    … …  
    74617439                        Element::GetInputListOnVertices(&GUrotm2[0],SealevelGUrotm2Enum);
    74627440                        Element::GetInputListOnVertices(&GUrotm3[0],SealevelGUrotm3Enum);
    7463                
     7441
    74647442                        for(int i=0;i<NUMVERTICES;i++){
    74657443                                if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    … …  
    74747452                                Element::GetInputListOnVertices(&GErotm2[0],SealevelGErotm2Enum);
    74757453                                Element::GetInputListOnVertices(&GErotm3[0],SealevelGErotm3Enum);
    7476                
     7454
    74777455                                for(int i=0;i<NUMVERTICES;i++){
    74787456                                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    … …  
    74957473        /*Create geoid: */
    74967474        for(int i=0;i<NUMVERTICES;i++)SealevelGrd[i]=UGrd[i]+RSLGrd[i];
    7497        
     7475
    74987476        /*Create inputs*/
    74997477        this->AddInput(SealevelGRDEnum,SealevelGrd,P1Enum);
    … …  
    75037481                this->AddInput(BedEastGRDEnum,EGrd,P1Enum);
    75047482        }
    7505 
    75067483
    75077484} /*}}}*/
    … …  
    75777554                }
    75787555        }
    7579        
     7556
    75807557        if(computeviscous){
    75817558                IssmDouble* grdfieldinterp=NULL;
    … …  
    76467623        }
    76477624
    7648 
    76497625} /*}}}*/
    76507626void       Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){/*{{{*/
    7651                
     7627
    76527628        IssmDouble S=0;
    76537629
    … …  
    76637639        longe=slgeom->longe[this->lid]/180*M_PI;
    76647640
    7665 
    76667641        /*recover total load: */
    76677642        if(loads->loads) S+=loads->loads[this->Sid()];
    76687643        if(loads->sealevelloads) S+=loads->sealevelloads[this->Sid()];
    7669        
     7644
    76707645        /* Perturbation terms for moment of inertia (moi_list):
    76717646         * computed analytically (see Wu & Peltier, eqs 10 & 32)
    … …  
    76797654}/*}}}*/
    76807655void       Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* grdloads,  SealevelGeometry* slgeom){/*{{{*/
    7681                
     7656
    76827657        IssmDouble  SA=0;
    76837658        IssmDouble* loads=NULL;
    … …  
    76917666        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    76927667        this->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
    7693        
     7668
    76947669        /*Initalize:*/
    76957670        for(int i=0;i<3;i++)dI_list[i]=0;
    … …  
    77237698}/*}}}*/
    77247699void       Tria::SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
    7725        
     7700
    77267701        if(slgeom->isoceanin[this->lid]){
    77277702                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    … …  
    77327707        }
    77337708
    7734 
    77357709} /*}}}*/
    77367710#endif
    … …  
    77417715        int interp;
    77427716        int type;
    7743        
     7717
    77447718        /*Branch according to whether we have a transient or not input: */
    77457719        type=this->inputs->GetInputObjectEnum(name);
    … …  
    77507724                this->InputServe(triainput);
    77517725                interp=triainput->GetInterpolation();
    7752                
     7726
    77537727                if (interp==P0Enum){
    77547728                        /*Update the value if this element belongs to the partition: */

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

    r26432 r26468  
    243243        /*If you reach this point, analysis has not been found*/
    244244        _error_("Could not find index of analysis " << EnumToStringx(analysis_enum) << " in list of FemModel analyses");
    245 
    246245
    247246}/*}}}*/
    … …  
    23552354        if(numonnodes) parameters->FindParam(&resultsonnodes,&numonnodes,SettingsResultsOnNodesEnum);
    23562355
    2357 
    23582356        /*Go through all requested output*/
    23592357        for(int i=0;i<numoutputs;i++){
    … …  
    23612359                output_enum   = StringToEnumx(output_string,false);
    23622360                isvec         = false;
    2363 
    23642361
    23652362                /*If string is not an enum, it is defined in output definitions*/
    … …  
    24812478                                        }
    24822479                                        break;
    2483 
    24842480
    24852481                                   /*Default is always Vector */
    … …  
    36583654        }
    36593655
    3660 
    36613656        /*Cleanup*/
    36623657        xDelete<IssmDouble>(P0inputs);
    … …  
    48124807        serial_active=active->ToMPISerial();
    48134808        delete active;
    4814 
    48154809
    48164810        /*Update node activation accordingly*/
    … …  
    54475441        bool refine;
    54485442
    5449 
    54505443        /*Fill variables*/
    54515444        switch(errorestimator_type){
    … …  
    55355528
    55365529        //itapopo esse metodo pode ser deletado
    5537 
    55385530
    55395531        /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/

  • TabularUnified issm/trunk-jpl/src/c/classes/GrdLoads.cpp

    r26227 r26468  
    1919GrdLoads::GrdLoads(int nel,SealevelGeometry* slgeom){ /*{{{*/
    2020
    21 
    2221        vloads=new Vector<IssmDouble>(nel);
    2322        for (int i=0;i<SLGEOM_NUMLOADS;i++) vsubloads[i]=new Vector<IssmDouble>(slgeom->nbar[i]);
    … …  
    2726
    2827        vsubsealevelloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
    29 
    3028
    3129}; /*}}}*/
    … …  
    5149                vsubloads[i]->Assemble();
    5250        }
    53        
     51
    5452        loads=vloads->ToMPISerial();
    5553        for (int i=0;i<SLGEOM_NUMLOADS;i++){
    … …  
    6563} /*}}}*/
    6664void GrdLoads::BroadcastSealevelLoads(void){ /*{{{*/
    67        
     65
    6866        sealevelloads=vsealevelloads->ToMPISerial();
    6967        subsealevelloads=vsubsealevelloads->ToMPISerial();

  • TabularUnified issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp

    r26160 r26468  
    223223        //}
    224224
    225 
    226225        //NEW??
    227226        //this->gradient->SetInput(interp,numindices,indices,values_in);

  • TabularUnified issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp

    r26323 r26468  
    390390/*}}}*/
    391391PentaInput* TransientInput::GetPentaInput(int offset){/*{{{*/
    392 
    393392
    394393        /*Check offset*/

  • TabularUnified issm/trunk-jpl/src/c/classes/IoModel.cpp

    r26432 r26468  
    612612                                record_name=xNew<char>(record_name_size+1);
    613613                                record_name[record_name_size]='\0';
    614                                
     614
    615615                                /*Read record_name: */
    616616                                if(fread(record_name,record_name_size*sizeof(char),1,fid)!=0){};
    617                                
     617
    618618                                _error_("error while looking in binary file. String \"" << record_name << "\" is a string of size "<<record_name_size);
    619619                        }
    … …  
    18381838                                matrix=array[i];
    18391839
    1840                                
    18411840                                //initialize transient input dataset:
    18421841                                TransientInput* transientinput=inputs->SetDatasetTransientInput(input_enum,i, times,N);
    … …  
    18501849                                        int        *vertexlids = xNew<int>(numvertices);
    18511850                                        int        *vertexsids = xNew<int>(numvertices);
    1852 
    18531851
    18541852                                        /*Recover vertices ids needed to initialize inputs*/
    … …  
    19181916                                        }
    19191917                                        else _error_("FetchDataToInput error message: row size of MatArray elements should be either numberofelements (+1) or numberofvertices (+1)");
    1920                                        
     1918
    19211919                                        xDelete<int>(vertexlids);
    19221920                                        xDelete<int>(vertexsids);

  • TabularUnified issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

    r26466 r26468  
    3434        element_in->FindParam(&this->law,FrictionLawEnum);
    3535        element_in->FindParam(&this->domaintype,DomainTypeEnum);
    36        
     36
    3737        /* Load VxBase and VyBase for this special case */
    3838        switch(this->domaintype){
    … …  
    6767        element_in->FindParam(&this->law,FrictionLawEnum);
    6868        element_in->FindParam(&this->domaintype,DomainTypeEnum);
    69        
     69
    7070        /* Load VxBase and VyBase for this special case */
    7171        switch(this->domaintype){
    … …  
    120120                        _error_("not supported");
    121121        }
    122 
    123122
    124123        /*Checks*/
    … …  
    278277        if(vmag==0. && (1./m-1.)<0.) alpha_complement=0.;
    279278        else alpha_complement= pow(vmag, 1.0/m-1.);
    280        
     279
    281280        /*Assign output pointers:*/
    282281        *palpha_complement=alpha_complement;
    … …  
    522521        Tpmp = element->TMeltingPoint(pressure);
    523522        deltaT = T-Tpmp;
    524 
    525523
    526524        /*Compute gamma*/
    … …  
    921919        if ((this->vz_input == NULL) || (this->apply_dim<3.)) vz = 0.0;
    922920        else this->vz_input->GetInputValue(&vz, gauss);
    923        
     921
    924922        if (this->apply_dim<2.) vy = 0.0;
    925923

  • TabularUnified issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r26294 r26468  
    752752        //IssmDouble eps_eff_averaged=0;
    753753        while(gauss_seg->next()){
    754                
     754
    755755                /*Compute zeta for gauss_seg point (0=surface, 1=base)*/
    756756                zeta=0.5*(gauss_seg->coord1+1);
    … …  
    773773                f[2]=(1-pow(zeta,n+1))*(1-pow(zeta,n+1));
    774774                f[3]=((n+1)/H)*pow(zeta,n) * ((n+1)/H)*pow(zeta,n);
    775        
     775
    776776                F[0]=H;
    777777                F[1]=H*(n+1)/(n+2);

  • TabularUnified issm/trunk-jpl/src/c/classes/Node.cpp

    r26245 r26468  
    654654                }
    655655                ug->SetValues(ssize,indices,values,INS_VAL);
    656                
     656
    657657                xDelete<IssmDouble>(values);
    658658                xDelete<int>(indices);
    … …  
    661661/*}}}*/
    662662void Node::VecReduce(Vector<IssmDouble>* uf, IssmDouble* local_ug,int* indices_ug){/*{{{*/
    663 
    664663
    665664        /*Only perform operation if not clone*/

  • TabularUnified issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp

    r26247 r26468  
    6969        int lower_row;
    7070
    71 
    7271        for (int i=0;i<SLGEOM_NUMLOADS;i++){
    7372                subelementmapping[i]=xNew<int>(localnel);
    … …  
    107106        /*Also, we'll need the barycentre associated areas:*/
    108107
    109 
    110108} /*}}}*/
    111109int SealevelGeometry::GEnum(int l){ /*{{{*/
    112        
     110
    113111        int output = -1;
    114112        switch(l){
    … …  
    122120} /*}}}*/
    123121int SealevelGeometry::GUEnum(int l){ /*{{{*/
    124        
     122
    125123        int output = -1;
    126124        switch(l){
    … …  
    134132} /*}}}*/
    135133int SealevelGeometry::GNEnum(int l){ /*{{{*/
    136        
     134
    137135        int output = -1;
    138136        switch(l){
    … …  
    146144} /*}}}*/
    147145int SealevelGeometry::GEEnum(int l){ /*{{{*/
    148        
     146
    149147        int output = -1;
    150148        switch(l){

  • TabularUnified issm/trunk-jpl/src/c/cores/WrapperPreCorePointerFromSolutionEnum.cpp

    r26222 r26468  
    2828                        break;
    2929        }
    30        
     30
    3131        /*Assign output pointer:*/
    3232        *psolutioncore=solutioncore;

  • TabularUnified issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp

    r23822 r26468  
    2828        if(VerboseSolution()) _printf0_("   computing velocities\n");
    2929        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
    30        
     30
    3131        bool is_schur_cg_solver = false;
    3232        #ifdef _HAVE_PETSC_

  • TabularUnified issm/trunk-jpl/src/c/cores/balancethickness_core.cpp

    r25680 r26468  
    1111
    1212void balancethickness_core(FemModel* femmodel){
    13 
    1413
    1514        /*recover parameters: */

  • TabularUnified issm/trunk-jpl/src/c/cores/control_core.cpp

    r26363 r26468  
    143143        /*Update control input*/
    144144        SetControlInputsFromVectorx(femmodel,X);
    145        
     145
    146146        /*solve forward: */
    147147        switch(solution_type){
    … …  
    152152                case StressbalanceSolutionEnum:{
    153153                        femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
    154                        
     154
    155155                        bool is_schur_cg_solver = false;
    156156                        #ifdef _HAVE_PETSC_
    … …  
    159159                        if(solver_type==FSSolverEnum) is_schur_cg_solver = true;
    160160                        #endif
    161                        
     161
    162162                        if(is_schur_cg_solver){
    163163                         solutionsequence_schurcg(femmodel);

  • TabularUnified issm/trunk-jpl/src/c/cores/dakota_core.cpp

    r24981 r26468  
    227227        if(VerboseQmu()) _printf0_("outputing results for this core:\n");
    228228        OutputResultsx(femmodel);
    229        
     229
    230230        /*Free ressources:*/
    231231        DakotaFree(&d_variables,&d_variables_descriptors,&responses_descriptors, d_numvariables, numresponsedescriptors);

  • TabularUnified issm/trunk-jpl/src/c/cores/esa_core.cpp

    r26055 r26468  
    114114                if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    115115        }
    116        
     116
    117117        /*End profiler*/
    118118        femmodel->profiler->Stop(ESACORE);

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

    r26121 r26468  
    4747                SolidEarthWaterUpdates(femmodel);
    4848
    49 
    5049        }
    5150        /*Using the Tws based Model*/
    5251        if (hydrology_model==HydrologyTwsEnum){
    5352                if(VerboseSolution()) _printf0_("   computing water column\n");
    54        
     53
    5554                femmodel->SetCurrentConfiguration(HydrologyTwsAnalysisEnum);
    56                
     55
    5756                /*save current tws  before updating:*/
    5857                InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
    59        
     58
    6059                /*grab tws from the hydrology.spcwatercolumn field input and update
    6160                 * the solution with it:*/
    … …  
    6867
    6968                delete ug;
    70        
     69
    7170        }
    7271
    … …  
    241240        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
    242241        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
    243        
     242
    244243        /*early return?:*/
    245244        if(!isgrd)return;

  • TabularUnified issm/trunk-jpl/src/c/cores/love_core.cpp

    r26272 r26468  
    2020                int nyi;
    2121                int starting_layer;
    22                
     22
    2323                LoveVariables(){  /*{{{*/
    2424                        g0=0;
    … …  
    245245        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
    246246        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
    247        
     247
    248248        g0=vars->g0;
    249249        r0=vars->r0;
    … …  
    392392                           fn=(deg*(deg+1.0));
    393393
    394 
    395394                           if(issolid==1){
    396395                           ny = 6;
    … …  
    652651                xmax=(matlitho->radius[i+1])/ra;
    653652
    654 
    655653                if (matlitho->issolid[i]){
    656654                        ny = 6;
    … …  
    662660                        one= -1.0;
    663661                }
    664 
    665662
    666663                for (int j = 0;j<ny;j++){
    … …  
    679676                        }
    680677                }
    681 
    682678
    683679                // Boundary Condition matrix - solid regions
    … …  
    712708        }
    713709
    714 
    715710        //-- Internal sphere: integration starts here rather than r=0 for numerical reasons
    716711        Innersphere_boundaryconditions<doubletype>(yi, starting_layer, deg, omega, femmodel, matlitho,vars);
    … …  
    856851
    857852                allgesv<doubletype>(&nyi, &nrhs, yilocal, &lda, ipiv, rhslocal, &ldb, &info);
    858                
     853
    859854                xDelete<int>(ipiv);
    860                
     855
    861856                if(VerboseModule() && info!=0){
    862857                        _printf0_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");
    … …  
    973968        IssmDouble* r          = matlitho->radius;
    974969        IssmDouble  r1,r2,ro, GG;
    975        
     970
    976971        /*outputs:*/
    977972        IssmDouble* EarthMass=NULL;
    … …  
    997992        nyi=6*(numlayers+1);
    998993        starting_layer=0;
    999        
     994
    1000995        return new LoveVariables(EarthMass,g0,r0,nyi,starting_layer);
    1001996
    … …  
    10191014        doubletype*  LoveHf = NULL;
    10201015        doubletype*  LoveLf = NULL;
    1021        
     1016
    10221017        doubletype*  LoveKernels= NULL;
    10231018        LoveVariables* vars=NULL;
    … …  
    10521047        femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum);
    10531048
    1054 
    10551049        /*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */
    10561050        LoveKf = xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
    … …  
    10981092        }
    10991093
    1100 
    11011094        xDelete<doubletype>(yi);
    11021095        xDelete<doubletype>(yi_righthandside);
    1103 
    11041096
    11051097        //Temporal love numbers
    … …  
    11451137                xDelete<IssmDouble>(LoveKtDouble);
    11461138                xDelete<IssmDouble>(LoveLtDouble);
    1147        
     1139
    11481140                /*Add into external results, no need to downcast, we can handle complexes/quad precision: */
    11491141                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0));
    … …  
    12261218        //              );
    12271219
    1228 
    12291220        /*Only when love_kernels is on*/
    12301221        //if (love_kernels==1) {
    … …  
    12321223                //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
    12331224        //}
    1234 
    12351225
    12361226        /*Add love matrices to results:*/
    … …  
    12411231        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0));
    12421232        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
    1243 
    1244 
    12451233
    12461234        //xDelete<IssmDouble>(LoveKr);

  • TabularUnified issm/trunk-jpl/src/c/cores/masstransport_core.cpp

    r26121 r26468  
    7777                }
    7878                SolidEarthIceUpdates(femmodel);
    79                
     79
    8080                femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
    8181                extrudefrombase_core(femmodel);
    … …  
    8484                femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
    8585                extrudefrombase_core(femmodel);
    86                
     86
    8787        }
    8888
    … …  
    112112        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
    113113        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
    114        
     114
    115115        /*early return?:*/
    116116        if(!isgrd)return;

  • TabularUnified issm/trunk-jpl/src/c/cores/movingfront_core.cpp

    r25883 r26468  
    1111
    1212void movingfront_core(FemModel* femmodel){
    13        
     13
    1414        /*Start profiler*/
    1515        femmodel->profiler->Start(MOVINGFRONTCORE);
    … …  
    4747                }
    4848        }
    49 
    5049
    5150        /* start the work from here */

  • TabularUnified issm/trunk-jpl/src/c/cores/oceantransport_core.cpp

    r26121 r26468  
    5151        /*profiler*/
    5252        femmodel->profiler->Stop(OCEANTRANSPORTCORE);
    53        
     53
    5454        /*free ressources:*/
    5555        delete ug;
    … …  
    6969        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
    7070        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
    71        
     71
    7272        /*early return?:*/
    7373        if(!isgrd)return;

  • TabularUnified issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

    r26296 r26468  
    88#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
    99#endif
    10 
    1110
    1211#include "./cores.h"
    … …  
    4443        /*Retrieve parameters:*/
    4544        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    46        
     45
    4746        /*Verbose: */
    4847        if(VerboseSolution()) _printf0_("   computing sea level change\n");
    … …  
    5352        /*run geometry core: */
    5453        slgeom=sealevelchange_geometry(femmodel);
    55        
     54
    5655        /*any external forcings?:*/
    5756        solidearthexternal_core(femmodel);
    … …  
    9796        int modelid=-1;
    9897        int  isexternal=0;
    99  
     98
    10099        /*parameters: */
    101100        IssmDouble          dt;
    … …  
    105104        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    106105        femmodel->parameters->FindParam(&isexternal,SolidearthIsExternalEnum);
    107        
     106
    108107        /*Early return:*/
    109108        if (!isexternal)return;
    … …  
    119118                GetVectorFromInputsx(&bedrocknorth,femmodel,BedNorthEnum,VertexSIdEnum);
    120119        }
    121        
     120
    122121        GetVectorFromInputsx(&geoid_rate,femmodel,SolidearthExternalGeoidRateEnum,VertexSIdEnum);
    123122        GetVectorFromInputsx(&bedrock_rate,femmodel,SolidearthExternalDisplacementUpRateEnum,VertexSIdEnum);
    … …  
    156155        /*Be very careful here, everything is well thought through, do not remove
    157156         * without taking big risks:*/
    158        
     157
    159158        /*parameters:*/
    160159        int  iscoupling;
    … …  
    199198        /*transer loads from basins to Earth for grd core computations. :*/
    200199        if(iscoupling){
    201        
     200
    202201                /*transfer ice thickness change load from basins to earth: */
    203202                TransferForcing(femmodel,DeltaIceThicknessEnum);
    … …  
    213212                }
    214213        }
    215        
     214
    216215}; /*}}}*/
    217216void              grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/
    … …  
    222221        GenericParam<BarystaticContributions*>* barycontribparam=NULL;
    223222        IssmDouble polarmotionvector[3]={0};
    224        
     223
    225224        GrdLoads*              loads=NULL;
    226225        Vector<IssmDouble>*    oldsealevelloads=NULL;
    … …  
    304303
    305304        if(VerboseSolution()) _printf0_("         starting  GRD convolutions\n");
    306        
     305
    307306        /*update viscous RSL:*/
    308307        for(Object* & object : femmodel->elements->objects){
    … …  
    341340                        element->SealevelchangeConvolution(sealevelpercpu, loads , polarmotionvector,slgeom);
    342341                }
    343                
     342
    344343                /*retrieve sea level average  and ocean area:*/
    345344                for(Object* & object : femmodel->elements->objects){
    … …  
    349348
    350349                loads->AssembleSealevelLoads();
    351        
     350
    352351                /*compute ocean areas:*/
    353352                if(!loads->sealevelloads){ //first time in the loop
    … …  
    357356                        if(scaleoceanarea) totaloceanarea=3.619e+14; // use true ocean area, m^2
    358357                }
    359        
     358
    360359                //Conserve ocean mass:
    361360                oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, totaloceanarea);
    … …  
    377376
    378377        deformation:
    379        
     378
    380379        if(VerboseSolution()) _printf0_("         deformation GRD convolutions\n");
    381380
    … …  
    403402                        oceanareas->Sum(&totaloceanarea); _assert_(totaloceanarea>0.);
    404403                        if(scaleoceanarea) totaloceanarea=3.619e+14; // use true ocean area, m^2
    405                                
     404
    406405                        //Conserve ocean mass
    407406                        //Note that here we create sea-level loads but they will not generate GRD as we have already run all the convolutions
    … …  
    454453        femmodel->parameters->FindParam(&isocean,TransientIsoceantransportEnum);
    455454        if(!isocean)return;
    456        
     455
    457456        /*Verbose: */
    458457        if(VerboseSolution()) _printf0_("         computing steric and dynamic sea level change\n");
    … …  
    477476        femmodel->parameters->SetParam(cumgmtslc,CumGmtslcEnum);
    478477        femmodel->parameters->SetParam(cumgmslc,CumGmslcEnum);
    479        
     478
    480479        /*Outputs some metrics:*/
    481480        femmodel->parameters->FindParam(&step,StepEnum);
    … …  
    496495/*}}}*/
    497496void              coupleroutput_core(FemModel* femmodel){  /*{{{*/
    498        
     497
    499498        /*parameters:*/
    500499        int iscoupling;
    … …  
    504503        femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
    505504        femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    506                
     505
    507506        if(iscoupling){
    508507                /*transfer sea level back to ice caps:*/
    … …  
    522521        IssmDouble          *xx     = NULL;
    523522        IssmDouble          *yy     = NULL;
    524        
     523
    525524        if(VerboseSolution()) _printf0_("         computing vertical deformation using Ivins model. \n");
    526525
    … …  
    541540        bedup = new Vector<IssmDouble>(gsize);
    542541        beduprate = new Vector<IssmDouble>(gsize);
    543        
     542
    544543        /*retrieve geometric information: */
    545544        VertexCoordinatesx(&xx,&yy,NULL,femmodel->vertices);
    … …  
    581580        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
    582581        nel=femmodel->elements->NumberOfElements();
    583                
     582
    584583        /*early return?:*/
    585584        if(grdmodel==IvinsEnum) return;
    … …  
    601600        femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel));
    602601        femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel));
    603 
    604602
    605603        #ifdef _ISSM_DEBUG_
    … …  
    612610        return;
    613611
    614 
    615612}/*}}}*/
    616613SealevelGeometry* sealevelchange_geometry(FemModel* femmodel) {  /*{{{*/
    … …  
    638635        femmodel->parameters->FindParam(&zze,&nel,ZzeEnum);
    639636        femmodel->parameters->FindParam(&areae,&nel,AreaeEnum);
    640        
     637
    641638        /*initialize SealevelloadMasks structure: */
    642639        slgeom=new SealevelGeometry(femmodel->elements->Size(),femmodel->vertices->Size());
    643        
     640
    644641        /*Verbose: */
    645642        if(VerboseSolution()) _printf0_("         computing sea level geometrical kernel and weight updates.\n");
    646        
    647        
     643
    648644        /*Run sealevel geometry routine for elements with full loading:*/
    649645        for(Object* & object : femmodel->elements->objects){
    … …  
    651647                element->SealevelchangeGeometryCentroidLoads(slgeom,xxe,yye,zze,areae);
    652648        }
    653        
     649
    654650        /*Initialize fractional loading mapping: */
    655651        slgeom->InitializeMappingsAndBarycentres();
    656652
    657        
    658653        /*Run sealevel geometry routine for elements with fractional loading:*/
    659654        for(Object* & object : femmodel->elements->objects){
    … …  
    671666        }
    672667
    673 
    674668        femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel));
    675669        femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel));
    … …  
    738732        vsealevelloadsvolume->PointwiseMult(loads->vsealevelloads,oceanareas);
    739733        vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas);
    740        
     734
    741735        vsealevelloadsvolume->Sum(&sealevelloadsaverage);
    742736        vsubsealevelloadsvolume->Sum(&subsealevelloadsaverage);
    743737        delete vsealevelloadsvolume;
    744738        delete vsubsealevelloadsvolume;
    745        
     739
    746740        return (sealevelloadsaverage+subsealevelloadsaverage)/totaloceanarea;
    747741} /*}}}*/
    … …  
    758752        IssmDouble      m1, m2, m3;
    759753        bool rotation=false;
    760        
     754
    761755        /*early return?:*/
    762756        femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    … …  
    786780        ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    787781        ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    788        
     782
    789783        ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    790784        ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    791        
     785
    792786        ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    793787        ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    … …  
    817811        int* indices=xNew<int>(nel);
    818812        for(int i=0;i<nel;i++)indices[i]=i;
    819        
     813
    820814        Vector<IssmDouble>* vloadcopy=new Vector<IssmDouble>(nel);
    821815        IssmDouble* loadcopy=xNew<IssmDouble>(nel);
    822        
     816
    823817        vloadcopy->SetValues(nel,indices,load,INS_VAL);
    824818        vloadcopy->Assemble();
    825 
    826819
    827820        if(subload){

  • TabularUnified issm/trunk-jpl/src/c/main/esmfbinders.cpp

    r26057 r26468  
    115115                                                Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
    116116                                                surface_input->GetInputAverage(&surface);
    117                        
     117
    118118                                                *(issmoutputs+f*numberofelements+i) = surface;
    119119

  • TabularUnified issm/trunk-jpl/src/c/main/issm.cpp

    r26047 r26468  
    1515        /*Initialize femmodel from arguments provided command line: */
    1616        FemModel *femmodel = new FemModel(argc,argv,comm_init);
    17        
     17
    1818        /*Need to know we are firing up from ISSM main, not a coupler driver like issm_slcp or issm_ocean:*/
    1919        femmodel->parameters->AddObject(new IntParam(IsSlcCouplingEnum,0));

  • TabularUnified issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp

    r23657 r26468  
    1212        int melt_parameterization;
    1313        femmodel->parameters->FindParam(&melt_parameterization,FrontalForcingsParamEnum);
    14        
     14
    1515        /*Calculate melting rate*/
    1616        switch(melt_parameterization){

  • TabularUnified issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp

    r26073 r26468  
    2727        }
    2828}
    29 
    30 
    3129
    3230void InputUpdateFromConstantx(FemModel* femmodel,IssmDouble constant, int name){
    … …  
    8482        else _error_("InputUpdateFromConstantx error message: type not supported yet!");
    8583
    86        
    8784}
    8885void InputUpdateFromConstantx(Inputs* inputs,Elements* elements,bool constant, int name){

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

    r26282 r26468  
    4949                descriptor=variables_descriptors[i];
    5050
    51 
    5251                /*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */
    5352                if (strncmp(descriptor,"scaled_",7)==0){
    … …  
    166165                        _printf0_("\n");
    167166                }
    168 
    169167
    170168                if (femmodel->inputs->GetInputObjectEnum(MasstransportSpcthicknessEnum)==DatasetInputEnum)
    … …  
    286284        }
    287285
    288 
    289286} /*}}}*/

  • TabularUnified issm/trunk-jpl/src/c/modules/MmeToInputFromIdx/MmeToInputFromIdx.cpp

    r26048 r26468  
    77#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
    88#endif
    9 
    109
    1110#include "../../shared/shared.h"

  • TabularUnified issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

    r25539 r26468  
    289289                        N_all[i] = N;
    290290
    291 
    292291                        for(Object* & object : elements->objects){
    293292                                Element* element=xDynamicCast<Element*>(object);

  • TabularUnified issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp

    r26252 r26468  
    717717                                xDelete<int>(cost_functions_weights_N);
    718718                                xDelete<IssmDouble*>(cost_functions_weights);
    719                        
     719
    720720                        /*}}}*/
    721721                        }

  • TabularUnified issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r26222 r26468  
    9090                parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum));
    9191                parameters->AddObject(iomodel->CopyConstantObject("md.transient.issampling",TransientIssamplingEnum));
    92                
     92
    9393                /*For stress balance only*/
    9494                parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum));

  • TabularUnified issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp

    r25627 r26468  
    2424        int         solutiontype;
    2525        char*       solutiontypestring      = NULL;
    26 
    2726
    2827        /*recover my_rank:*/
    … …  
    7170                                char* root=NULL;
    7271                                char* modelname=NULL;
    73                                
     72
    7473                                femmodel->parameters->FindParam(&currEvalId,QmuCurrEvalIdEnum);
    7574                                femmodel->parameters->FindParam(&statistics,QmuStatisticsEnum);

  • TabularUnified issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp

    r25711 r26468  
    2020        IssmDouble* dmatfield = NULL;
    2121        int*        imatfield = NULL;
    22                        
     22
    2323        //type of the returned field:
    2424        int type;
    … …  
    109109        int range,lower_row,upper_row;
    110110        int nfilesperdirectory;
    111        
     111
    112112        /*intermediary:*/
    113113        IssmDouble* doublemat=NULL;
    … …  
    237237                        char* field=fields[f];
    238238                        meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
    239                        
     239
    240240                        if(meanxtype[f]==1){
    241241                                meanxsize[f]=1;
    … …  
    246246                                        timemean+=scalar/nsteps;
    247247                                }
    248                                        
     248
    249249                                /*Figure out max and min of time means: */
    250250                                if(i==(lower_row+1)){
    … …  
    334334                                IssmDouble*  allmax=xNew<IssmDouble>(size);
    335335                                IssmDouble*  allmin=xNew<IssmDouble>(size);
    336                                
     336
    337337                                ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    338338                                ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    … …  
    344344                } /*}}}*/
    345345        }
    346        
     346
    347347        /*Now do the same for the time mean fields:*/
    348348        for (int f=0;f<nfields;f++){
    … …  
    377377                        ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    378378
    379                        
    380379                        /*Store for later use in histogram computation:*/
    381380                        maxmeans[f]=allmax;
    … …  
    488487                        char* field=fields[f];
    489488                        meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
    490                        
     489
    491490                        if(meanxtype[f]==1){
    492491                                IssmDouble timemean=0;
    … …  
    496495                                        timemean+=scalar/nsteps;
    497496                                }
    498                                        
     497
    499498                                /*Figure out max and min of time means: */
    500499                                if(i==(lower_row+1)){
    … …  
    531530                                        IssmDouble* ma=maxmeans[f];
    532531                                        IssmDouble* mi=minmeans[f];
    533                                        
     532
    534533                                        for (int k=0;k<doublematsize;k++){
    535534                                                IssmDouble scalar=timemean[k];
    … …  
    541540                                }
    542541                                else{
    543                                        
     542
    544543                                        IssmDouble* localhistogram=timehistogram[f];
    545544                                        IssmDouble* ma=maxmeans[f];
    546545                                        IssmDouble* mi=minmeans[f];
    547                                        
     546
    548547                                        for (int k=0;k<doublematsize;k++){
    549548                                                IssmDouble scalar=timemean[k];
    … …  
    600599                                IssmDouble* histo=histogram[counter];
    601600                                IssmDouble* allhisto=xNew<IssmDouble>(size*nbins);
    602                                
     601
    603602                                ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    604603                                xDelete<IssmDouble>(histo);
    … …  
    620619        }
    621620        _printf0_("Start aggregating time mean histogram:\n");
    622        
     621
    623622        /*Now do the same for the time mean fields:*/
    624623        for (int f=0;f<nfields;f++){
    … …  
    684683        int range,lower_row,upper_row;
    685684        int nfilesperdirectory;
    686        
     685
    687686        /*intermediary:*/
    688687        IssmDouble* doublemat=NULL;
    … …  
    806805                        }
    807806                }
    808                
     807
    809808                /*Deal with time mean: */
    810809                for (int f=0;f<nfields;f++){
    … …  
    917916                                IssmDouble*  sumx=xNew<IssmDouble>(size);
    918917                                IssmDouble*  sumx2=xNew<IssmDouble>(size);
    919                                
     918
    920919                                ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    921920                                ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    … …  
    980979                        IssmDouble*  sumx=xNew<IssmDouble>(size);
    981980                        IssmDouble*  sumx2=xNew<IssmDouble>(size);
    982                                
     981
    983982                        ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    984983                        ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    … …  
    10211020        int* indices=NULL;
    10221021        int  nindices;
    1023        
     1022
    10241023        /*intermediary:*/
    10251024        IssmDouble* doublemat=NULL;
    … …  
    11431142                                if(my_rank==0){
    11441143                                        char fieldname[1000];
    1145                                        
     1144
    11461145                                        sprintf(fieldname,"%s%s",fields[f],"Samples");
    11471146                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
    … …  
    11521151                                x=xs[counter];
    11531152                                allx=xNew<IssmDouble>(nsamples*nindices);
    1154                                
     1153
    11551154                                MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
    1156                                
     1155
    11571156                                /*add to results:*/
    11581157                                if(my_rank==0){
    … …  
    11831182
    11841183        FemModel* femmodel=new FemModel();
    1185        
     1184
    11861185        /*Some parameters that will allow us to use the OutputResultsx module:*/
    11871186        parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
    … …  
    13271326                results   = new Results();
    13281327                parameters=new Parameters();
    1329                
     1328
    13301329                //solution type:
    13311330                parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum));
    1332        
     1331
    13331332                //root  directory
    13341333                directory=xNew<char>(strlen(argv[2])+1);
    … …  
    13591358                ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm);
    13601359
    1361 
    13621360                iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics");
    13631361                for (int i=1;i<=numstatistics;i++){
    … …  
    13671365                        int   nsamples;
    13681366                        _printf0_("Dealing with qmu statistical computation #" << i << "\n");
    1369                
     1367
    13701368                        sprintf(string,"md.qmu.statistics.method(%i).name",i);
    13711369                        iomodel->FindConstant(&name,string);
    … …  
    13911389                                iomodel->FetchData(&indices,&dummy,&nindices,string);
    13921390                                parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices));
    1393                
     1391
    13941392                                ComputeSampleSeries(parameters,results,color,statcomm);
    13951393                        }
    … …  
    14101408        /*output results:*/
    14111409        OutputStatistics(parameters,results,color,statcomm);
    1412        
     1410
    14131411        /*all meet here: */
    14141412        ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n");

  • TabularUnified issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

    r26271 r26468  
    16181618
    16191619        //// MELT, PERCOLATION AND REFREEZE
    1620        
     1620
    16211621        // initialize refreeze, runoff, flxDn and dW vectors [kg]
    16221622        IssmDouble* F = xNewZeroInit<IssmDouble>(n);
    … …  
    17071707                                for(int l=i;(l<n && d[l]>=dPHC-Dtol);l++) depthice=depthice+dz[l];
    17081708                        }
    1709                        
     1709
    17101710                        // break loop if there is no meltwater and if depth is > mw_depth
    17111711                        if (fabs(inM) < Wtol && i > X){

  • TabularUnified issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp

    r23400 r26468  
    3232  IssmDouble Pmax = 0.6;
    3333  IssmDouble inv_twelve=1./12.;
    34  
     34
    3535  sconv=(rho_water/rho_ice);            //rhow_rain/rhoi
    3636
    … …  
    9191         IssmDouble coeff5=3.8e-05;             // polynomial
    9292         IssmDouble coeff6=6.0e-07;             // fit
    93    
     93
    9494         if(tstar>=temp_rain)
    9595          frac_solid = 0.0;
    … …  
    100100                                         +tstar*(coeff3+tstar*(coeff4+tstar*(coeff5+tstar*coeff6))));
    101101         }
    102          
     102
    103103         snowfall=snowfall+precip_star*frac_solid*inv_twelve;
    104104  }
    … …  
    108108  if(snowfall<0.0) snowfall=0.0;   // correction of
    109109  if(rainfall<0.0) rainfall=0.0;   // negative values
    110    
     110
    111111  if(rainfall<=(Pmax*snowfall)){
    112112          if((rainfall+beta1*pdd)<=(Pmax*snowfall)) {
    … …  
    126126          runoff     = smelt+rainfall-Pmax*snowfall;
    127127  }
    128    
     128
    129129  saccu = precip;       
    130130

  • TabularUnified issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp

    r25549 r26468  
    7676                Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
    7777                femmodel->profiler->Stop(SOLVER);
    78        
     78
    7979                Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
    8080
    … …  
    9898                        }
    9999                }
    100                
     100
    101101                /*Increase count: */
    102102                count++;

  • TabularUnified issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp

    r25710 r26468  
    3333        IssmPDouble                             t1,t2;                                  /* Time measurement for bottleneck analysis */
    3434
    35        
    3635        double tmp1,tmp2,tmp3;
    3736        int tmpi;
    … …  
    3938
    4039        int noIt;
    41 
    42 
    4340
    4441        int precond = 0;
    … …  
    4946        PetscBool flag,flg;
    5047        #endif
    51 
    5248
    5349        char ksp_type[50];
    … …  
    7167        #endif
    7268
    73        
    7469        if(precond){
    7570                _printf0_("Running WITH preconditioner\n");
    … …  
    7873        }
    7974
    80 
    8175        /*Initialize output*/
    8276        Vector<IssmDouble>* out_uf=new Vector<IssmDouble>(uf0);
    83        
     77
    8478        /* Extract block matrices from the saddle point matrix */
    8579        /* [ A   B ] = Kff
    … …  
    9791        MatGetSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT);
    9892        #endif
    99        
     93
    10094        /* Extract preconditioner matrix on the pressure space*/
    10195        #if PETSC_VERSION_GT(3,8,0)
    … …  
    10599        #endif
    106100
    107 
    108101        /* Get number of velocity / pressure nodes */
    109102        MatGetSize(B,&nu,&np);
    … …  
    117110        VecGetSubVector(out_uf->pvector->vector,isv,&uold);
    118111        VecGetSubVector(out_uf->pvector->vector,isp,&p);
    119 
    120112
    121113        /* Set up intermediaries */
    … …  
    133125        VecDuplicate(p,&wnew);VecSet(wnew,0.0);
    134126        VecDuplicate(uold,&chi);VecSet(chi,0.0);
    135        
     127
    136128        /* Get global RHS (for each block sub-problem respectively)*/
    137129        VecGetSubVector(pf,isv,&f1);
    … …  
    143135        /* a(u0,v) = f1(v)-b(p0,v)  i.e.  Au0 = F1-Bp0 */
    144136        /* u0 = u_DIR on \Gamma_DIR */
    145        
     137
    146138        /* Create KSP context */
    147139        KSPCreate(IssmComm::GetComm(),&kspu);
    … …  
    174166                _error_("Suggested KSP method not implemented yet!\n");
    175167        }
    176        
     168
    177169        KSPSetInitialGuessNonzero(kspu,PETSC_TRUE);
    178        
     170
    179171        /*Strong rel. residual tolerance needed when solving for the velocity update.
    180172         * This is because ISSM uses the dimensional equations, so the initial guess chi = 0
    181173         * results in a very high initial residual (~ 1e19)*/
    182174        KSPSetTolerances(kspu,ELLTOL,PETSC_DEFAULT,PETSC_DEFAULT,maxiter); //maxiter
    183 
    184 
    185175
    186176        KSPGetPC(kspu,&pcu);
    … …  
    210200        KSPSetUp(kspu);
    211201
    212        
    213202        /* Create RHS */
    214203        /* RHS = F1-B * pold */
    … …  
    227216                VecNorm(resu,NORM_2,&tmp5);
    228217
    229 
    230         }
    231 
     218        }
    232219
    233220        /* Go solve Au0 = F1-Bp0*/
    … …  
    235222        KSPGetIterationNumber(kspu,&noIt);
    236223
    237 
    238 
    239224        if (VerboseConvergence())
    240225        {
    241        
     226
    242227        KSPGetIterationNumber(kspu,&tmpi);
    243228        VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.);
    … …  
    249234        }
    250235
    251 
    252 
    253236        /* Set up u_new */
    254237        VecDuplicate(uold,&unew);VecCopy(uold,unew);
    255238        VecAssemblyBegin(unew);VecAssemblyEnd(unew);
    256239
    257 
    258 
    259240        /* ------------------------------------------------------------- */
    260241
    261242        /*Get initial residual*/
    262243        /*(1/mu(x) * g0, q) = b(q,u0) - (f2,q)  i.e.  IP * g0 = BT * u0 - F2*/
    263        
     244
    264245        /* Create KSP context */
    265246        KSPCreate(IssmComm::GetComm(),&kspip);
    … …  
    269250        KSPSetOperators(kspip,IP,IP,DIFFERENT_NONZERO_PATTERN);
    270251        #endif
    271        
     252
    272253        /* Create RHS */
    273254        /* RHS = BT * uold - F2 */
    274255        VecScale(uold,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(uold,-1.);
    275256
    276 
    277257        /* Set KSP & PC options */
    278258        KSPSetType(kspip,KSPGMRES);
    279259        KSPSetInitialGuessNonzero(kspip,PETSC_TRUE);
    280        
    281        
     260
    282261        KSPGetPC(kspip,&pcp);
    283262        PCSetType(pcp,PCJACOBI);
    … …  
    294273        }
    295274
    296 
    297 
    298275        /* Go solve */
    299276        KSPSolve(kspip,rhsp,gold);
    300277
    301 
    302 
    303278        if (VerboseConvergence())
    304279        {
    305        
     280
    306281        KSPGetResidualNorm(kspip,&tmp1);
    307282        VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.);
    … …  
    311286        _printf0_("||IP*g00 - rhsp||_euc: " << tmp4 <<"\n||Jac(-1)*(IP*g00-rhsp)||_euc: " << tmp5 << "\n||IP*gn0-rhsp||_euc: " << tmp2<< "\n||Jac(-1)*(IP*gn0-rhsp)||_euc: " << tmp1 << "\nIteration number: "<<tmpi<<"\n");   
    312287        }
    313 
    314288
    315289        /*Initial residual*/
    … …  
    318292        _printf0_("inner product norm g0: "<<initRnorm<<"\n");
    319293
    320 
    321294        /*Iterate only if inital residual large enough*/
    322295        if(initRnorm > 1e-5)
    … …  
    326299        VecDuplicate(gold,&gnew);VecCopy(gold,gnew);
    327300        VecAssemblyBegin(gnew);VecAssemblyEnd(gnew);
    328 
    329301
    330302        /* ------------------------------------------------------------ */
    … …  
    334306        VecDuplicate(gold,&wold);VecCopy(gold,wold);
    335307        VecAssemblyBegin(wold);VecAssemblyEnd(wold);
    336 
    337308
    338309        /* Count number of iterations */
    … …  
    351322                VecScale(wold,1.);MatMult(B,wold,rhsu);VecScale(wold,1.);
    352323                VecSet(chi,0.0);
    353        
    354                
     324
    355325                if(VerboseConvergence())
    356326                {
    357327                VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.);
    358328                VecNorm(tmpu,NORM_2,&tmp4);
    359                        
     329
    360330                KSPInitialResidual(kspu,chi,tmpu,tmpu2,resu,rhsu);
    361331                VecNorm(resu,NORM_2,&tmp5);
    … …  
    363333                }
    364334
    365                
    366335                        KSPSolve(kspu,rhsu,chi);
    367                
    368        
    369                
    370                
     336
    371337                if (VerboseConvergence())
    372338                {
    373339                VecNorm(chi,NORM_2,&tmp1);
    374340                KSPGetResidualNorm(kspu,&tmp2);
    375        
     341
    376342                VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.);
    377343                VecNorm(tmpu,NORM_2,&tmp3);
    378                
    379                
     344
    380345                KSPGetIterationNumber(kspu,&tmpi);
    381346                _printf0_("||chi_nk||_euc: "<< tmp1 << "\n||A*chi_0k - rhsu||_euc: "<<tmp4<< "\n||P(-1)*(A*chi_0k-rhsu)||_euc: " << tmp5 << "\n||A*chi_nk - rhsu||_euc: "<< tmp3 <<"\n||P(-1)*(A*chi_nk - rhsu)||_euc: " << tmp2 <<"\nIteration Number: " << tmpi <<"\n");
    382347                }
    383348
    384 
    385349                /* ---------------------------------------------------------- */
    386 
    387350
    388351                /*Set step size*/
    … …  
    390353                MatMult(IP,gold,tmpp);
    391354                VecDot(tmpp,wold,&rho);
    392                
     355
    393356                MatMult(BT,chi,tmpp);
    394357                VecDot(tmpp,wold,&tmpScalar);
    395358                rho = rho/tmpScalar;
    396359
    397 
    398360                /* ---------------------------------------------------------- */
    399 
    400361
    401362                /*Pressure update*/
    402363                /*p(m+1) = pm + rhom * wm*/
    403364                VecAXPY(p,-1.*rho,wold);
    404 
    405365
    406366                /*Velocity update*/
    … …  
    408368                VecWAXPY(unew,rho,chi,uold);
    409369                VecNorm(unew,NORM_2,&tmpScalar);
    410 
    411370
    412371                if (VerboseConvergence())
    … …  
    421380                _printf0_("Incompressibility norm: " << tmp7 <<"\n");
    422381                }
    423 
    424 
    425382
    426383                /* ---------------------------------------------------------- */
    … …  
    432389                VecWAXPY(rhsp,-1.*rho,tmpp,tmpp2);
    433390                KSPSolve(kspip,rhsp,gnew);
    434                
    435 
    436391
    437392                /* ---------------------------------------------------------- */
    438393
    439394                MatMult(IP,gnew,tmpp);
    440                
     395
    441396                VecDot(tmpp,gnew,&gamma);
    442397                rnorm = sqrt(gamma);
    … …  
    454409                        _printf0_("L2-Norm g: "<< rnorm << "\n");
    455410                }
    456        
     411
    457412                /*Break prematurely if solver doesn't reach desired tolerance in reasonable time frame*/
    458413                if(count > 10./TOL)   
    … …  
    461416                 break;
    462417                }
    463        
    464418
    465419                /* ---------------------------------------------------------- */
    466 
    467420
    468421                /*Directional update*/
    … …  
    475428                /* Assign new to old iterates */
    476429                VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold);
    477                
     430
    478431                count++;
    479432        }
    … …  
    487440        *puf=out_uf;
    488441
    489 
    490442        /* Cleanup */
    491443        KSPDestroy(&kspu);KSPDestroy(&kspip);
    492444
    493445        MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);MatDestroy(&IP);
    494        
     446
    495447        VecDestroy(&p);VecDestroy(&uold);VecDestroy(&unew);VecDestroy(&rhsu);VecDestroy(&rhsp);
    496448        VecDestroy(&gold);VecDestroy(&gnew);VecDestroy(&wold);VecDestroy(&wnew);VecDestroy(&chi);
    … …  
    522474        double rnorm1, rnorm2;
    523475
    524 
    525476        if(VerboseModule()) _printf0_("   checking convergence\n");
    526477
    … …  
    548499        MatGetSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT);
    549500        #endif
    550        
     501
    551502                /*no. of free nodes in velocity/pressure space*/
    552503                MatGetSize(B,&n_u,&n_p);
    … …  
    560511                VecGetSubVector(uf->pvector->vector,isv,&u);
    561512                VecGetSubVector(uf->pvector->vector,isp,&p);
    562                
    563513
    564514                /*Extract values of the RHS corresponding to the first/second block*/
    … …  
    573523                VecDuplicate(p,&res2);VecSet(res2,1.0);
    574524
    575 
    576525        /*Display solver caracteristics*/
    577526        if (VerboseConvergence()){
    578                
     527
    579528                /*Calculate res1 = A*u + B*p - f1*/
    580529                VecScale(f1,-1.);MatMultAdd(A,u,f1,tmp);MatMultAdd(B,p,tmp,res1);VecScale(f1,-1.);
    581530                /*Calculate res2 = B^T * u - f2*/
    582531                VecScale(f2,-1.);MatMultAdd(BT,u,f2,res2);VecScale(f2,-1.);
    583 
    584532
    585533                /*compute norm(res1), norm(res2), norm(F) and residue*/
    … …  
    597545        VecRestoreSubVector(uf->pvector->vector,isv,&u);
    598546        VecRestoreSubVector(uf->pvector->vector,isp,&p);
    599        
     547
    600548        /*Extract values corresponding to velocity/pressure on the old solution*/
    601549        VecGetSubVector(old_uf->pvector->vector,isv,&uold);
    602550        VecGetSubVector(old_uf->pvector->vector,isp,&pold);
    603                
    604551
    605552        /*Force equilibrium (Mandatory)*/
    … …  
    609556        /*Calculate res2 = B^T * uold - f2*/
    610557        VecScale(f2,-1.);MatMultAdd(BT,uold,f2,res2);VecScale(f2,-1.);
    611        
     558
    612559        /*compute norm(res1), norm(res2), norm(F) and residue*/
    613560        VecNorm(res1,NORM_2,&rnorm1);VecNorm(res2,NORM_2,&rnorm2);
    … …  
    626573        VecRestoreSubVector(old_uf->pvector->vector,isv,&uold);
    627574        VecRestoreSubVector(old_uf->pvector->vector,isp,&pold);
    628        
    629 
    630575
    631576        //print
    … …  
    721666        int  count=0;
    722667        bool converged=false;
    723        
     668
    724669        /*Start non-linear iteration using input velocity: */
    725670        GetSolutionFromInputsx(&ug,femmodel);
    … …  
    729674        InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
    730675        InputUpdateFromSolutionx(femmodel,ug);
    731        
     676
    732677        for(;;){
    733678
    … …  
    747692                PetscOptionsGetInt(NULL,PETSC_NULL,"-schur_pc",&precond,NULL);
    748693                #endif
    749        
     694
    750695                StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); 
    751        
     696
    752697                /* Construct Schur preconditioner matrix or mass matrix depending on input */
    753698                if(precond)
    … …  
    760705                        }
    761706                }else{
    762                        
     707
    763708                        for(Object* & object : femmodel->elements->objects){
    764709                                Element* element2=xDynamicCast<Element*>(object);
    … …  
    769714                }
    770715
    771                        
    772716                delete analysis;
    773        
     717
    774718                /*Obtain index sets for velocity and pressure components */
    775719                IS isv = NULL;
    … …  
    785729                Kff->Assemble();
    786730
    787        
    788731                /*Solve*/
    789732                femmodel->profiler->Start(SOLVER);
    790733                _assert_(Kff->type==PetscMatType);
    791 
    792734
    793735                SchurCGSolver(&uf,

  • TabularUnified issm/trunk-jpl/src/c/toolkits/codipack/ampi_interface.cpp

    r23350 r26468  
    1212#error "Cannot compile without MeDiPack and AdjointMpi"
    1313#endif
    14 

  • TabularUnified issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp

    r26433 r26468  
    318318  xDelete(valueX);
    319319  xDelete(indexX);
    320  
     320
    321321  delete data;
    322322}
Note: See TracChangeset for help on using the changeset viewer.