Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 26468) @@ -150,5 +150,5 @@ }/*}}}*/ ElementMatrix* AdjointHorizAnalysis::CreateKMatrixHO(Element* element){/*{{{*/ - + /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 26468) @@ -70,5 +70,4 @@ element->SetElementInput(inputs,DamageFEnum,0.); } - /*What input do I need to run my damage evolution model?*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 26468) @@ -535,5 +535,4 @@ int iseplthickcomp; - /*Skip if we don't want to compute thicknesses*/ femmodel->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum); @@ -558,5 +557,4 @@ femmodel->parameters->FindParam(&gravity,ConstantsGEnum); femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - for(Object* & object : femmodel->elements->objects){ Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 26468) @@ -562,5 +562,5 @@ /*Skip if water or ice shelf element*/ if(element->IsAllFloating()) return; - + /*Intermediary*/ IssmDouble phi_0, phi_m, p_i; Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 26468) @@ -307,5 +307,5 @@ /*Add input to the element: */ element->AddInput(WatercolumnEnum,watercolumn,element->GetElementType()); - + /*Also update the hydrological loads for the sealevel core: */ IssmDouble* oldwatercolumn = xNew(numnodes); Index: /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp (revision 26468) @@ -42,5 +42,4 @@ iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",WatercolumnEnum); - }/*}}}*/ void HydrologyTwsAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -85,5 +84,5 @@ int numnodes = element->GetNumberOfNodes(); IssmDouble* watercolumn = xNew(numnodes); - + /*Use the dof list to index into the solution vector: */ for(int i=0;iparameters->FindParam(TimesteppingTimeStepEnum); - for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 26468) @@ -234,5 +234,4 @@ if(grdmodel==IvinsEnum) inputs->SetTransientInput(TransientAccumulatedDeltaIceThicknessEnum,NULL,0); - }/*}}}*/ void MasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -816,5 +815,5 @@ element->FindParam(&count,SealevelchangeRunCountEnum); } - + /*Fetch dof list and allocate solution vector*/ int *doflist = NULL; @@ -912,5 +911,4 @@ element->AddBasalInput(SurfaceEnum,newsurface,P1Enum); element->AddBasalInput(BaseEnum,newbase,P1Enum); - /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/analyses/OceantransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/OceantransportAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/OceantransportAnalysis.cpp (revision 26468) @@ -62,5 +62,4 @@ iomodel->FetchDataToInput(inputs,elements,"md.initialization.str",StrEnum); - }/*}}}*/ void OceantransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -75,5 +74,5 @@ IssmDouble modelid; int nummodels; - + /*create double param, not int param, because Dakota will be updating it as * a double potentially: */ Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26468) @@ -49,12 +49,11 @@ iomodel->FetchDataToInput(inputs,elements,"md.solidearth.external.geoid",SolidearthExternalGeoidRateEnum); - /*Resolve Mmes using the modelid, if necessary:*/ if (inputs->GetInputObjectEnum(SolidearthExternalDisplacementEastRateEnum)==DatasetInputEnum){ int modelid; - + /*retrieve model id: */ iomodel->FetchData(&modelid,"md.solidearth.external.modelid"); - + /*replace dataset of forcings with only one, the modelid'th:*/ MmeToInputFromIdx(inputs,elements,modelid,SolidearthExternalDisplacementNorthRateEnum, P1Enum); @@ -74,5 +73,4 @@ iomodel->ConstantToInput(inputs,elements,0.,DeltaIceThicknessEnum,P1Enum); iomodel->ConstantToInput(inputs,elements,0.,DeltaBottomPressureEnum,P1Enum); - }/*}}}*/ @@ -182,5 +180,5 @@ xDelete(partitionice); } - + iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro"); parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro)); @@ -208,5 +206,5 @@ BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel); parameters->AddObject(new GenericParam(barystaticcontributions,BarystaticContributionsEnum)); - + /*Deal with external multi-model ensembles: {{{*/ if(isexternal){ @@ -425,6 +423,4 @@ xDelete(displs); - - /*Avoid singularity at 0: */ G_gravi[0]=G_gravi[1]; @@ -436,6 +432,4 @@ } } - - /*Reinterpolate viscoelastic green kernels onto a regular gridded time @@ -484,5 +478,5 @@ } } - + } else if(elastic){ Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 26468) @@ -232,5 +232,4 @@ parameters->AddObject(new IntParam(SmbStepsPerStepEnum,smbslices)); - parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum)); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26468) @@ -1109,5 +1109,4 @@ #endif - if(is_schur_cg_solver) solutionsequence_schurcg(femmodel); @@ -2096,5 +2095,5 @@ /*Use the dof list to index into the solution vector: */ for(i=0;iTransformSolutionCoord(&values[0],XYEnum); @@ -2134,5 +2133,5 @@ } element->AddBasalInput(VelEnum,vel,element->GetElementType()); - + /*Free ressources:*/ xDelete(vel); @@ -2754,7 +2753,7 @@ /*MLHO*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/ - + //_error_("Mono Layer Higher-Order called, not fully tested. If you are sure of using it, comment this line."); - + /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; @@ -2818,5 +2817,5 @@ gauss = element->NewGauss(2); } - + /* Start looping on the number of gaussian points: */ while(gauss->next()){ @@ -2850,5 +2849,4 @@ xDelete(basis); return Ke; - //itapopo OLD - testing above @@ -2873,5 +2871,5 @@ //ElementMatrix* Ke = basalelement->NewElementMatrix(MLHOApproximationEnum); ElementMatrix* KeSSA = CreateKMatrixSSAFriction(basalelement); //only to get K11 and K33 - + /*Fetch number of nodes and dof for this finite element*/ //int numnodes = basalelement->GetNumberOfNodes(); @@ -2969,5 +2967,5 @@ 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] )*(n+1)/(n+2);//K14 - + Ke->values[(4*i+1)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity[1]*thickness*( 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] @@ -2985,5 +2983,5 @@ 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] )*2*pow(n+1,2)/((2*n+3)*(n+2));//K24 - + Ke->values[(4*i+2)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity[0]*thickness*( 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i] @@ -3017,5 +3015,5 @@ } } - + /*Transform Coordinate System*/ //basalelement->TransformStiffnessMatrixCoord(Ke,XYMLHOEnum); @@ -3051,5 +3049,5 @@ default: _error_("mesh "<GetInputDerivativeValue(&slope[0],xyz_list,gauss); n_input->GetInputValue(&n,gauss); - + for(int i=0;ivalues[i*4+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F1 @@ -3150,5 +3148,5 @@ element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); element->NodalFunctions(basis,gauss); - + b = base-sealevel; // ice base shifted by the sea level s = thickness+b; // ice surface shifted by the sea level @@ -3249,5 +3247,5 @@ IssmDouble* vel = xNew(numnodes); IssmDouble* n = xNew(numnodes); - + /*Use the dof list to index into the solution vector: */ for(i=0;iAddBasalInput(VxShearEnum,vshx,element->GetElementType()); element->AddBasalInput(VyShearEnum,vshy,element->GetElementType()); - + /*Add vx and vy as inputs to the tria element (base velocities): */ element->AddBasalInput(VxBaseEnum,vbx,element->GetElementType()); element->AddBasalInput(VyBaseEnum,vby,element->GetElementType()); - + /*Add vx and vy as inputs to the tria element (surface velocities): */ element->AddBasalInput(VxSurfaceEnum,vsx,element->GetElementType()); @@ -3298,5 +3296,5 @@ vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2); } - + /*Get Vz and compute vel (vertically averaged vel)*/ basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); @@ -3372,5 +3370,5 @@ solution->SetValues(numdof,doflist,values,INS_VAL); - + /*Free ressources:*/ delete gauss; @@ -4406,5 +4404,4 @@ if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} - /* Start looping on the number of gaussian points: */ Gauss* gauss = element->NewGauss(5); @@ -4414,5 +4411,4 @@ element->NodalFunctionsPressure(pbasis,gauss); element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - if(dim==3 || dim==2){ @@ -5485,5 +5481,5 @@ element->NodalFunctionsPressure(pbasis,gauss); element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - + for(int i=0;iFetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum); #endif - /*Add basal forcings to compute melt rate*/ @@ -583,5 +582,4 @@ } - #ifdef INWOOVZ IssmDouble* thickness = xNew(numnodes); @@ -624,5 +622,4 @@ xDelete(thickness); #endif - /*Now Compute vel*/ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 26468) @@ -86,5 +86,5 @@ /*Methods*/ void AmrBamg::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/ - + /*Delete previous mesh and keep the entire mesh*/ if(this->elementslist) xDelete(this->elementslist); @@ -99,5 +99,5 @@ }/*}}}*/ void AmrBamg::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/ - + /*Get the entire mesh*/ *elementslist_out = this->elementslist; @@ -128,5 +128,5 @@ }/*}}}*/ void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int** pdatalist,IssmDouble** pxylist,int** pelementslist){/*{{{*/ - + /*Intermediaries*/ BamgGeom* geomout=new BamgGeom(); @@ -142,5 +142,5 @@ this->options->field = field; this->options->hmaxVertices = hmaxVertices; - + /*remesh*/ if(this->previousmesh){ Index: /issm/trunk-jpl/src/c/classes/AmrNeopz.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrNeopz.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/AmrNeopz.cpp (revision 26468) @@ -137,5 +137,5 @@ /*Mesh refinement methods*/ void AmrNeopz::SetMesh(int** elementslist_in,IssmDouble** x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/*{{{*/ - + /*Delete previous mesh and keep the entire mesh*/ if(this->elementslist) xDelete(this->elementslist); @@ -150,5 +150,5 @@ }/*}}}*/ void AmrNeopz::GetMesh(int** elementslist_out,IssmDouble** x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/*{{{*/ - + /*Get the entire mesh*/ *elementslist_out = this->elementslist; @@ -394,5 +394,5 @@ if(type>1) break; } - + /*Verify and return*/ if(type>1) type=2; @@ -416,5 +416,5 @@ count=1; - + while(count>0){ count=0; @@ -438,5 +438,5 @@ } } - + /*refine the element if requested*/ if(type==2){gmesh->Element(i)->Divide(sons); count++;} @@ -926,5 +926,5 @@ this->fathermesh = dynamic_cast(TPZPersistenceManager::ReadFromFile()); TPZPersistenceManager::CloseRead(); - + TPZPersistenceManager::OpenRead(previousmeshfile); this->previousmesh = dynamic_cast(TPZPersistenceManager::ReadFromFile()); @@ -933,9 +933,9 @@ if(!amrifstream.is_open()) _error_("amr ifstream is not open!"); amrifstream.seekg(0); - + this->sid2index.clear(); this->index2sid.clear(); this->specialelementsindex.clear(); - + amrifstream>>size; for(int i=0;isid2index.push_back(value); } - + amrifstream>>size; for(int i=0;iindex2sid.push_back(value); } - + amrifstream>>size; for(int i=0;iindex2sid.size()>0){ amrofstream << this->index2sid.size() << std::endl; Index: /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp (revision 26468) @@ -66,5 +66,5 @@ IssmDouble sumice,sumhydro,sumocean; - + ice->Assemble(); hydro->Assemble(); @@ -90,5 +90,4 @@ cumocean->Sum(&sumocean); - return sumice+sumhydro+sumocean; @@ -100,8 +99,7 @@ cumhydro->AXPY(hydro,1); - } /*}}}*/ void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/ - + int id; if(nice){ @@ -128,14 +126,12 @@ ocean->SetValue(0,oceanvalue,ADD_VAL); } - } /*}}}*/ void BarystaticContributions::Reset(){ /*{{{*/ - ice->Set(0.); hydro->Set(0.); ocean->Set(0.); - + } /*}}}*/ void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp (revision 26468) @@ -96,5 +96,4 @@ IssmDouble Cfsurfacelogvel::Response(FemModel* femmodel){/*{{{*/ - /*recover time parameters: */ IssmDouble time; Index: /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 26468) @@ -124,5 +124,5 @@ if(VerboseQmu()) _printf0_("compute dakota responses:\n"); femmodel->DakotaResponsesx(responses,responses_descriptors,numresponsedescriptors,numFns); - + /*Output results for this iteration: */ if(VerboseQmu()) _printf0_("output results for this iteration: \n"); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26468) @@ -1286,5 +1286,4 @@ /*}}}*/ void Element::GetVectorFromInputs(Vector* vector,int input_enum,int type){/*{{{*/ - switch(type){ @@ -1680,8 +1679,8 @@ if(M==1)times[0]=0; if(M==2)for(int t=0;tSetTransientInput(vector_enum,times,N); TransientInput* transientinput = inputs->GetTransientInput(vector_enum); - + for(int t=0;tGetNumberOfVertices(); - - this->GetVerticesSidList(&sidlist[0]); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26468) @@ -1940,5 +1940,4 @@ } - /*Some checks in debugging mode*/ _assert_(s1>=0 && s1<=1.); @@ -2055,5 +2054,4 @@ _error_("case not possible"); } - /*Some checks in debugging mode*/ @@ -2377,5 +2375,4 @@ if(start==+1 && !IsOnSurface()) return; - /*Get original input*/ Input* input = this->GetInput(enum_type); @@ -2738,5 +2735,5 @@ IssmDouble movingfrontvy[NUMVERTICES]; IssmDouble vel; - + /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -4033,5 +4030,4 @@ } - /*Some checks in debugging mode*/ _assert_(s1>=0 && s1<=1.); @@ -4147,5 +4143,4 @@ _error_("case not possible"); } - /*Some checks in debugging mode*/ @@ -4899,5 +4894,4 @@ else if(type==TransientInputEnum){ - IssmDouble* steps=NULL; int nsteps; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26468) @@ -323,5 +323,4 @@ IssmDouble epse_2,groundedice,bed,sealevel; // added sealevel - /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -337,6 +336,4 @@ Input* n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input* sl_input = this->GetInput(SealevelEnum); _assert_(sl_input); - - /* Start looping on the number of vertices: */ @@ -413,5 +410,4 @@ int crevasse_opening_stress; - /*reset if no ice in element*/ if(!this->IsIceInElement()){ @@ -426,5 +422,4 @@ return; } - /*retrieve the type of crevasse_opening_stress*/ @@ -1321,5 +1316,5 @@ area=this->GetAreaSpherical(); vareae->SetValue(this->sid,area,INS_VAL); - + /*in addition, put in in the inputs:*/ this->inputs->SetDoubleInput(AreaEnum,this->lid,area); @@ -1350,5 +1345,5 @@ vlate->SetValue(this->sid,late,INS_VAL); vareae->SetValue(this->sid,area,INS_VAL); - + return; } @@ -1766,5 +1761,5 @@ /*}}}*/ void Tria::GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){/*{{{*/ - + /*Computeportion of the element that is grounded*/ bool trapezeisnegative=true; //default value @@ -1777,5 +1772,5 @@ IssmDouble loadweights_g[NUMVERTICES]; IssmDouble total_weight=0; - + _assert_(!xIsNan(gl[0])); _assert_(!xIsNan(gl[1])); @@ -1820,9 +1815,8 @@ if(trapezeisnegative) phi=1-f1*f2; else phi=f1*f2; - /*Compute weights:*/ gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2); - + total_weight=0; for(int i=0;iparameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); @@ -1940,5 +1934,5 @@ Element::GetInputListOnVertices(&levelset[0],levelsetenum); if(flip1)for(int i=0;iGetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset); @@ -1951,5 +1945,4 @@ } /*}}}*/ void Tria::GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum){ /*{{{*/ - bool istrapeze1, istrapeze2; @@ -2019,5 +2012,4 @@ } - //If everyone is negative, no need to calculate any fraction if (levelset1[0]<=0 && levelset1[1]<=0 && levelset1[2]<=0 && levelset2[0]<=0 && levelset2[1]<=0 && levelset2[2]<=0) { @@ -2028,8 +2020,7 @@ return; } - + /*recover planet radius:*/ this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); - //If just one levelset is all negative, just take the partitioning of the other, no interaction between them @@ -2047,8 +2038,6 @@ } - this->GetFractionGeometry(&weights1[0],&phi1,&point1,&d,&e,&istrapeze1,levelset1); this->GetFractionGeometry(&weights2[0],&phi2,&point2,&f,&g,&istrapeze2,levelset2); - //Early return if levelsets are not independent @@ -2069,5 +2058,4 @@ } - ::GetVerticesCoordinates(&xyz0[0][0],vertices,NUMVERTICES); // initial triangle @@ -2117,5 +2105,4 @@ h2= (d*(f+e-1))/(g*(e-d) +f*d); } - //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]) @@ -2227,5 +2214,5 @@ } /*}}}*/ } - + }/*}}}*/ else{ /*{{{*/ @@ -2498,5 +2485,5 @@ }/*}}}*/ Input* Tria::GetInput(int inputenum){/*{{{*/ - + /*Get Input from dataset*/ if(this->iscollapsed){ @@ -4081,5 +4068,5 @@ IssmDouble movingfrontvy[NUMVERTICES]; IssmDouble vel; - + /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -4332,5 +4319,5 @@ } for(i=0;iparameters->FindParam(&gsize,MeshNumberofverticesEnum); this->parameters->FindParam(&yts,ConstantsYtsEnum); - + /*recover gia solution parameters: */ int cross_section_shape; @@ -6158,5 +6145,5 @@ /*recover material parameters: */ IssmDouble rho_ice = FindParam(MaterialsRhoIceEnum); - + /*recover mantle and lithosphere material properties:*/ int numlayers=litho->numlayers; @@ -6179,5 +6166,4 @@ int numtimes; this->GetInputAveragesUpToCurrentTime(TransientAccumulatedDeltaIceThicknessEnum,&hes,×,&numtimes,currenttime); - /*pull area of this Tria: */ @@ -6274,5 +6260,5 @@ IssmDouble* H_viscoelastic_precomputed=NULL; #endif - + /*viscoelastic green function:*/ int index; @@ -6342,5 +6328,5 @@ parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL); - + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL); @@ -6373,5 +6359,4 @@ } - constant=3/rho_earth/planetarea; @@ -6393,5 +6378,5 @@ IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0))); IssmDouble longe= atan2(yye[e],xxe[e]); - + lati=latitude[i]; longi=longitude[i]; @@ -6404,5 +6389,4 @@ _assert_(index>=0 && indexparameters->FindParam(&computeice,TransientIsmasstransportEnum); @@ -6671,5 +6654,5 @@ slgeom->nsubel[SLGEOM_OCEAN]++; } - + /*early return if we are not on an ice sheet , and we are not requesting *hydrology or bottom pressure loads :*/ @@ -6748,5 +6731,5 @@ } } - + /*Deal with ice loads if we are on grounded ice:*/ if(isice && !isoceanonly && computeice){ @@ -6791,5 +6774,5 @@ } } - + } /*}}}*/ @@ -6798,5 +6781,5 @@ /* Classic buildup of load weights, centroids and areas *for elements which are fully inside a mask. * At the same time, we'll tag the elements that are fractionally only inside a mask*/ - + IssmDouble loadweights[3]={0}; IssmDouble area,loadarea; @@ -6808,5 +6791,5 @@ IssmDouble constant; IssmDouble nanconstant=NAN; - + /*get vertex and area information:*/ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -6821,8 +6804,8 @@ this->AddInput(SealevelBarystaticOceanLongbarEnum,&longbar,P0Enum); #endif - + if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; - + this->GetNodalWeightsAndAreaAndCentroidsFromLeveset(&loadweightsocean[0],&loadareaocean,&latbar, &longbar, slgeom->late[this->lid], slgeom->longe[this->lid], area, MaskOceanLevelsetEnum); slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=loadareaocean; @@ -6885,5 +6868,5 @@ this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum); this->AddInput(SealevelBarystaticHydroAreaEnum,&loadarea,P0Enum); - + this->AddInput(SealevelBarystaticHydroLatbarEnum,&latbar,P0Enum); this->AddInput(SealevelBarystaticHydroLongbarEnum,&longbar,P0Enum); @@ -6891,5 +6874,5 @@ #endif } - + } /*}}}*/ @@ -6927,5 +6910,5 @@ IssmDouble** GEsubel=NULL; #endif - + /*viscoelastic green function:*/ int index; @@ -7010,5 +6993,5 @@ } } - + //Allocate: for(int l=0;lparameters->FindParam(&viscous,SolidearthSettingsViscousEnum); - + if(viscous){ this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - + this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum); @@ -7165,5 +7148,4 @@ this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&dummy); } - bool foundtime=false; @@ -7194,5 +7176,4 @@ viscousindex=newindex+offset; - this->parameters->SetParam(viscousindex,SealevelchangeViscousIndexEnum); this->parameters->SetParam(viscoustimes,viscousnumsteps,SealevelchangeViscousTimesEnum); @@ -7201,5 +7182,4 @@ xDelete(viscoustimes); } - } @@ -7254,9 +7234,9 @@ bslchydro = -slgeom->LoadArea[SLGEOM_WATER][this->lid]*Wavg; bslcbp = -slgeom->LoadArea[SLGEOM_OCEAN][this->lid]*BPavg; - + _assert_(!xIsNan(bslcice)); _assert_(!xIsNan(bslchydro)); _assert_(!xIsNan(bslcbp)); - + /*Plug values into subelement load vector:*/ if(slgeom->issubelement[SLGEOM_ICE][this->lid]){ @@ -7289,5 +7269,5 @@ IssmDouble* Gsub[SLGEOM_NUMLOADS]; bool computefuture=false; - + bool sal = false; bool viscous = false; @@ -7299,5 +7279,5 @@ IssmDouble Grotm2[3]; IssmDouble Grotm3[3]; - + this->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum); this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum); @@ -7318,5 +7298,5 @@ Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum); Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); - + for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ @@ -7333,5 +7313,5 @@ IssmDouble oceanarea=0; IssmDouble rho_water; - + this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); @@ -7352,5 +7332,4 @@ } else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL); - /*add ocean area into a global oceanareas vector:*/ @@ -7362,5 +7341,5 @@ } } - + return; } /*}}}*/ @@ -7409,6 +7388,5 @@ this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); this->parameters->FindParam(&planethasocean,SolidearthSettingsGrdOceanEnum); - - + if(sal){ @@ -7450,5 +7428,5 @@ Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum); Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); - + for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ @@ -7461,5 +7439,5 @@ Element::GetInputListOnVertices(&GUrotm2[0],SealevelGUrotm2Enum); Element::GetInputListOnVertices(&GUrotm3[0],SealevelGUrotm3Enum); - + for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ @@ -7474,5 +7452,5 @@ Element::GetInputListOnVertices(&GErotm2[0],SealevelGErotm2Enum); Element::GetInputListOnVertices(&GErotm3[0],SealevelGErotm3Enum); - + for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ @@ -7495,5 +7473,5 @@ /*Create geoid: */ for(int i=0;iAddInput(SealevelGRDEnum,SealevelGrd,P1Enum); @@ -7503,5 +7481,4 @@ this->AddInput(BedEastGRDEnum,EGrd,P1Enum); } - } /*}}}*/ @@ -7577,5 +7554,5 @@ } } - + if(computeviscous){ IssmDouble* grdfieldinterp=NULL; @@ -7646,8 +7623,7 @@ } - } /*}}}*/ void Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){/*{{{*/ - + IssmDouble S=0; @@ -7663,9 +7639,8 @@ longe=slgeom->longe[this->lid]/180*M_PI; - /*recover total load: */ if(loads->loads) S+=loads->loads[this->Sid()]; if(loads->sealevelloads) S+=loads->sealevelloads[this->Sid()]; - + /* Perturbation terms for moment of inertia (moi_list): * computed analytically (see Wu & Peltier, eqs 10 & 32) @@ -7679,5 +7654,5 @@ }/*}}}*/ void Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* grdloads, SealevelGeometry* slgeom){/*{{{*/ - + IssmDouble SA=0; IssmDouble* loads=NULL; @@ -7691,5 +7666,5 @@ this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); this->parameters->FindParam(&re,SolidearthPlanetRadiusEnum); - + /*Initalize:*/ for(int i=0;i<3;i++)dI_list[i]=0; @@ -7723,5 +7698,5 @@ }/*}}}*/ void Tria::SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/ - + if(slgeom->isoceanin[this->lid]){ if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ @@ -7732,5 +7707,4 @@ } - } /*}}}*/ #endif @@ -7741,5 +7715,5 @@ int interp; int type; - + /*Branch according to whether we have a transient or not input: */ type=this->inputs->GetInputObjectEnum(name); @@ -7750,5 +7724,5 @@ this->InputServe(triainput); interp=triainput->GetInterpolation(); - + if (interp==P0Enum){ /*Update the value if this element belongs to the partition: */ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26468) @@ -243,5 +243,4 @@ /*If you reach this point, analysis has not been found*/ _error_("Could not find index of analysis " << EnumToStringx(analysis_enum) << " in list of FemModel analyses"); - }/*}}}*/ @@ -2355,5 +2354,4 @@ if(numonnodes) parameters->FindParam(&resultsonnodes,&numonnodes,SettingsResultsOnNodesEnum); - /*Go through all requested output*/ for(int i=0;i(P0inputs); @@ -4812,5 +4807,4 @@ serial_active=active->ToMPISerial(); delete active; - /*Update node activation accordingly*/ @@ -5447,5 +5441,4 @@ bool refine; - /*Fill variables*/ switch(errorestimator_type){ @@ -5535,5 +5528,4 @@ //itapopo esse metodo pode ser deletado - /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/ Index: /issm/trunk-jpl/src/c/classes/GrdLoads.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/GrdLoads.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/GrdLoads.cpp (revision 26468) @@ -19,5 +19,4 @@ GrdLoads::GrdLoads(int nel,SealevelGeometry* slgeom){ /*{{{*/ - vloads=new Vector(nel); for (int i=0;i(slgeom->nbar[i]); @@ -27,5 +26,4 @@ vsubsealevelloads=new Vector(slgeom->nbar[SLGEOM_OCEAN]); - }; /*}}}*/ @@ -51,5 +49,5 @@ vsubloads[i]->Assemble(); } - + loads=vloads->ToMPISerial(); for (int i=0;iToMPISerial(); subsealevelloads=vsubsealevelloads->ToMPISerial(); Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 26468) @@ -223,5 +223,4 @@ //} - //NEW?? //this->gradient->SetInput(interp,numindices,indices,values_in); Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp (revision 26468) @@ -390,5 +390,4 @@ /*}}}*/ PentaInput* TransientInput::GetPentaInput(int offset){/*{{{*/ - /*Check offset*/ Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 26468) @@ -612,8 +612,8 @@ record_name=xNew(record_name_size+1); record_name[record_name_size]='\0'; - + /*Read record_name: */ if(fread(record_name,record_name_size*sizeof(char),1,fid)!=0){}; - + _error_("error while looking in binary file. String \"" << record_name << "\" is a string of size "<SetDatasetTransientInput(input_enum,i, times,N); @@ -1850,5 +1849,4 @@ int *vertexlids = xNew(numvertices); int *vertexsids = xNew(numvertices); - /*Recover vertices ids needed to initialize inputs*/ @@ -1918,5 +1916,5 @@ } else _error_("FetchDataToInput error message: row size of MatArray elements should be either numberofelements (+1) or numberofvertices (+1)"); - + xDelete(vertexlids); xDelete(vertexsids); Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26468) @@ -34,5 +34,5 @@ element_in->FindParam(&this->law,FrictionLawEnum); element_in->FindParam(&this->domaintype,DomainTypeEnum); - + /* Load VxBase and VyBase for this special case */ switch(this->domaintype){ @@ -67,5 +67,5 @@ element_in->FindParam(&this->law,FrictionLawEnum); element_in->FindParam(&this->domaintype,DomainTypeEnum); - + /* Load VxBase and VyBase for this special case */ switch(this->domaintype){ @@ -120,5 +120,4 @@ _error_("not supported"); } - /*Checks*/ @@ -278,5 +277,5 @@ if(vmag==0. && (1./m-1.)<0.) alpha_complement=0.; else alpha_complement= pow(vmag, 1.0/m-1.); - + /*Assign output pointers:*/ *palpha_complement=alpha_complement; @@ -522,5 +521,4 @@ Tpmp = element->TMeltingPoint(pressure); deltaT = T-Tpmp; - /*Compute gamma*/ @@ -921,5 +919,5 @@ if ((this->vz_input == NULL) || (this->apply_dim<3.)) vz = 0.0; else this->vz_input->GetInputValue(&vz, gauss); - + if (this->apply_dim<2.) vy = 0.0; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 26468) @@ -752,5 +752,5 @@ //IssmDouble eps_eff_averaged=0; while(gauss_seg->next()){ - + /*Compute zeta for gauss_seg point (0=surface, 1=base)*/ zeta=0.5*(gauss_seg->coord1+1); @@ -773,5 +773,5 @@ f[2]=(1-pow(zeta,n+1))*(1-pow(zeta,n+1)); f[3]=((n+1)/H)*pow(zeta,n) * ((n+1)/H)*pow(zeta,n); - + F[0]=H; F[1]=H*(n+1)/(n+2); Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 26468) @@ -654,5 +654,5 @@ } ug->SetValues(ssize,indices,values,INS_VAL); - + xDelete(values); xDelete(indices); @@ -661,5 +661,4 @@ /*}}}*/ void Node::VecReduce(Vector* uf, IssmDouble* local_ug,int* indices_ug){/*{{{*/ - /*Only perform operation if not clone*/ Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp (revision 26468) @@ -69,5 +69,4 @@ int lower_row; - for (int i=0;i(localnel); @@ -107,8 +106,7 @@ /*Also, we'll need the barycentre associated areas:*/ - } /*}}}*/ int SealevelGeometry::GEnum(int l){ /*{{{*/ - + int output = -1; switch(l){ @@ -122,5 +120,5 @@ } /*}}}*/ int SealevelGeometry::GUEnum(int l){ /*{{{*/ - + int output = -1; switch(l){ @@ -134,5 +132,5 @@ } /*}}}*/ int SealevelGeometry::GNEnum(int l){ /*{{{*/ - + int output = -1; switch(l){ @@ -146,5 +144,5 @@ } /*}}}*/ int SealevelGeometry::GEEnum(int l){ /*{{{*/ - + int output = -1; switch(l){ Index: /issm/trunk-jpl/src/c/cores/WrapperPreCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/WrapperPreCorePointerFromSolutionEnum.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/WrapperPreCorePointerFromSolutionEnum.cpp (revision 26468) @@ -28,5 +28,5 @@ break; } - + /*Assign output pointer:*/ *psolutioncore=solutioncore; Index: /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp (revision 26468) @@ -28,5 +28,5 @@ if(VerboseSolution()) _printf0_(" computing velocities\n"); femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); - + bool is_schur_cg_solver = false; #ifdef _HAVE_PETSC_ Index: /issm/trunk-jpl/src/c/cores/balancethickness_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancethickness_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/balancethickness_core.cpp (revision 26468) @@ -11,5 +11,4 @@ void balancethickness_core(FemModel* femmodel){ - /*recover parameters: */ Index: /issm/trunk-jpl/src/c/cores/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 26468) @@ -143,5 +143,5 @@ /*Update control input*/ SetControlInputsFromVectorx(femmodel,X); - + /*solve forward: */ switch(solution_type){ @@ -152,5 +152,5 @@ case StressbalanceSolutionEnum:{ femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); - + bool is_schur_cg_solver = false; #ifdef _HAVE_PETSC_ @@ -159,5 +159,5 @@ if(solver_type==FSSolverEnum) is_schur_cg_solver = true; #endif - + if(is_schur_cg_solver){ solutionsequence_schurcg(femmodel); Index: /issm/trunk-jpl/src/c/cores/dakota_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/dakota_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/dakota_core.cpp (revision 26468) @@ -227,5 +227,5 @@ if(VerboseQmu()) _printf0_("outputing results for this core:\n"); OutputResultsx(femmodel); - + /*Free ressources:*/ DakotaFree(&d_variables,&d_variables_descriptors,&responses_descriptors, d_numvariables, numresponsedescriptors); Index: /issm/trunk-jpl/src/c/cores/esa_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 26468) @@ -114,5 +114,5 @@ if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} } - + /*End profiler*/ femmodel->profiler->Stop(ESACORE); Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 26468) @@ -47,15 +47,14 @@ SolidEarthWaterUpdates(femmodel); - } /*Using the Tws based Model*/ if (hydrology_model==HydrologyTwsEnum){ if(VerboseSolution()) _printf0_(" computing water column\n"); - + femmodel->SetCurrentConfiguration(HydrologyTwsAnalysisEnum); - + /*save current tws before updating:*/ InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum); - + /*grab tws from the hydrology.spcwatercolumn field input and update * the solution with it:*/ @@ -68,5 +67,5 @@ delete ug; - + } @@ -241,5 +240,5 @@ femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); - + /*early return?:*/ if(!isgrd)return; Index: /issm/trunk-jpl/src/c/cores/love_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 26468) @@ -20,5 +20,5 @@ int nyi; int starting_layer; - + LoveVariables(){ /*{{{*/ g0=0; @@ -245,5 +245,5 @@ femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); - + g0=vars->g0; r0=vars->r0; @@ -392,5 +392,4 @@ fn=(deg*(deg+1.0)); - if(issolid==1){ ny = 6; @@ -652,5 +651,4 @@ xmax=(matlitho->radius[i+1])/ra; - if (matlitho->issolid[i]){ ny = 6; @@ -662,5 +660,4 @@ one= -1.0; } - for (int j = 0;j(yi, starting_layer, deg, omega, femmodel, matlitho,vars); @@ -856,7 +851,7 @@ allgesv(&nyi, &nrhs, yilocal, &lda, ipiv, rhslocal, &ldb, &info); - + xDelete(ipiv); - + if(VerboseModule() && info!=0){ _printf0_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system"); @@ -973,5 +968,5 @@ IssmDouble* r = matlitho->radius; IssmDouble r1,r2,ro, GG; - + /*outputs:*/ IssmDouble* EarthMass=NULL; @@ -997,5 +992,5 @@ nyi=6*(numlayers+1); starting_layer=0; - + return new LoveVariables(EarthMass,g0,r0,nyi,starting_layer); @@ -1019,5 +1014,5 @@ doubletype* LoveHf = NULL; doubletype* LoveLf = NULL; - + doubletype* LoveKernels= NULL; LoveVariables* vars=NULL; @@ -1052,5 +1047,4 @@ femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum); - /*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */ LoveKf = xNewZeroInit(nfreq*(sh_nmax+1)); @@ -1098,8 +1092,6 @@ } - xDelete(yi); xDelete(yi_righthandside); - //Temporal love numbers @@ -1145,5 +1137,5 @@ xDelete(LoveKtDouble); xDelete(LoveLtDouble); - + /*Add into external results, no need to downcast, we can handle complexes/quad precision: */ femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0)); @@ -1226,5 +1218,4 @@ // ); - /*Only when love_kernels is on*/ //if (love_kernels==1) { @@ -1232,5 +1223,4 @@ //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0)); //} - /*Add love matrices to results:*/ @@ -1241,6 +1231,4 @@ //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0)); //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0)); - - //xDelete(LoveKr); Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 26468) @@ -77,5 +77,5 @@ } SolidEarthIceUpdates(femmodel); - + femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); @@ -84,5 +84,5 @@ femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); - + } @@ -112,5 +112,5 @@ femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); - + /*early return?:*/ if(!isgrd)return; Index: /issm/trunk-jpl/src/c/cores/movingfront_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/movingfront_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/movingfront_core.cpp (revision 26468) @@ -11,5 +11,5 @@ void movingfront_core(FemModel* femmodel){ - + /*Start profiler*/ femmodel->profiler->Start(MOVINGFRONTCORE); @@ -47,5 +47,4 @@ } } - /* start the work from here */ Index: /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp (revision 26468) @@ -51,5 +51,5 @@ /*profiler*/ femmodel->profiler->Stop(OCEANTRANSPORTCORE); - + /*free ressources:*/ delete ug; @@ -69,5 +69,5 @@ femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); - + /*early return?:*/ if(!isgrd)return; Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26468) @@ -8,5 +8,4 @@ #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif - #include "./cores.h" @@ -44,5 +43,5 @@ /*Retrieve parameters:*/ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - + /*Verbose: */ if(VerboseSolution()) _printf0_(" computing sea level change\n"); @@ -53,5 +52,5 @@ /*run geometry core: */ slgeom=sealevelchange_geometry(femmodel); - + /*any external forcings?:*/ solidearthexternal_core(femmodel); @@ -97,5 +96,5 @@ int modelid=-1; int isexternal=0; - + /*parameters: */ IssmDouble dt; @@ -105,5 +104,5 @@ femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); femmodel->parameters->FindParam(&isexternal,SolidearthIsExternalEnum); - + /*Early return:*/ if (!isexternal)return; @@ -119,5 +118,5 @@ GetVectorFromInputsx(&bedrocknorth,femmodel,BedNorthEnum,VertexSIdEnum); } - + GetVectorFromInputsx(&geoid_rate,femmodel,SolidearthExternalGeoidRateEnum,VertexSIdEnum); GetVectorFromInputsx(&bedrock_rate,femmodel,SolidearthExternalDisplacementUpRateEnum,VertexSIdEnum); @@ -156,5 +155,5 @@ /*Be very careful here, everything is well thought through, do not remove * without taking big risks:*/ - + /*parameters:*/ int iscoupling; @@ -199,5 +198,5 @@ /*transer loads from basins to Earth for grd core computations. :*/ if(iscoupling){ - + /*transfer ice thickness change load from basins to earth: */ TransferForcing(femmodel,DeltaIceThicknessEnum); @@ -213,5 +212,5 @@ } } - + }; /*}}}*/ void grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/ @@ -222,5 +221,5 @@ GenericParam* barycontribparam=NULL; IssmDouble polarmotionvector[3]={0}; - + GrdLoads* loads=NULL; Vector* oldsealevelloads=NULL; @@ -304,5 +303,5 @@ if(VerboseSolution()) _printf0_(" starting GRD convolutions\n"); - + /*update viscous RSL:*/ for(Object* & object : femmodel->elements->objects){ @@ -341,5 +340,5 @@ element->SealevelchangeConvolution(sealevelpercpu, loads , polarmotionvector,slgeom); } - + /*retrieve sea level average and ocean area:*/ for(Object* & object : femmodel->elements->objects){ @@ -349,5 +348,5 @@ loads->AssembleSealevelLoads(); - + /*compute ocean areas:*/ if(!loads->sealevelloads){ //first time in the loop @@ -357,5 +356,5 @@ if(scaleoceanarea) totaloceanarea=3.619e+14; // use true ocean area, m^2 } - + //Conserve ocean mass: oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, totaloceanarea); @@ -377,5 +376,5 @@ deformation: - + if(VerboseSolution()) _printf0_(" deformation GRD convolutions\n"); @@ -403,5 +402,5 @@ oceanareas->Sum(&totaloceanarea); _assert_(totaloceanarea>0.); if(scaleoceanarea) totaloceanarea=3.619e+14; // use true ocean area, m^2 - + //Conserve ocean mass //Note that here we create sea-level loads but they will not generate GRD as we have already run all the convolutions @@ -454,5 +453,5 @@ femmodel->parameters->FindParam(&isocean,TransientIsoceantransportEnum); if(!isocean)return; - + /*Verbose: */ if(VerboseSolution()) _printf0_(" computing steric and dynamic sea level change\n"); @@ -477,5 +476,5 @@ femmodel->parameters->SetParam(cumgmtslc,CumGmtslcEnum); femmodel->parameters->SetParam(cumgmslc,CumGmslcEnum); - + /*Outputs some metrics:*/ femmodel->parameters->FindParam(&step,StepEnum); @@ -496,5 +495,5 @@ /*}}}*/ void coupleroutput_core(FemModel* femmodel){ /*{{{*/ - + /*parameters:*/ int iscoupling; @@ -504,5 +503,5 @@ femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum); femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - + if(iscoupling){ /*transfer sea level back to ice caps:*/ @@ -522,5 +521,5 @@ IssmDouble *xx = NULL; IssmDouble *yy = NULL; - + if(VerboseSolution()) _printf0_(" computing vertical deformation using Ivins model. \n"); @@ -541,5 +540,5 @@ bedup = new Vector(gsize); beduprate = new Vector(gsize); - + /*retrieve geometric information: */ VertexCoordinatesx(&xx,&yy,NULL,femmodel->vertices); @@ -581,5 +580,5 @@ femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); nel=femmodel->elements->NumberOfElements(); - + /*early return?:*/ if(grdmodel==IvinsEnum) return; @@ -601,5 +600,4 @@ femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel)); femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel)); - #ifdef _ISSM_DEBUG_ @@ -612,5 +610,4 @@ return; - }/*}}}*/ SealevelGeometry* sealevelchange_geometry(FemModel* femmodel) { /*{{{*/ @@ -638,12 +635,11 @@ femmodel->parameters->FindParam(&zze,&nel,ZzeEnum); femmodel->parameters->FindParam(&areae,&nel,AreaeEnum); - + /*initialize SealevelloadMasks structure: */ slgeom=new SealevelGeometry(femmodel->elements->Size(),femmodel->vertices->Size()); - + /*Verbose: */ if(VerboseSolution()) _printf0_(" computing sea level geometrical kernel and weight updates.\n"); - - + /*Run sealevel geometry routine for elements with full loading:*/ for(Object* & object : femmodel->elements->objects){ @@ -651,9 +647,8 @@ element->SealevelchangeGeometryCentroidLoads(slgeom,xxe,yye,zze,areae); } - + /*Initialize fractional loading mapping: */ slgeom->InitializeMappingsAndBarycentres(); - /*Run sealevel geometry routine for elements with fractional loading:*/ for(Object* & object : femmodel->elements->objects){ @@ -671,5 +666,4 @@ } - femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel)); femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel)); @@ -738,10 +732,10 @@ vsealevelloadsvolume->PointwiseMult(loads->vsealevelloads,oceanareas); vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas); - + vsealevelloadsvolume->Sum(&sealevelloadsaverage); vsubsealevelloadsvolume->Sum(&subsealevelloadsaverage); delete vsealevelloadsvolume; delete vsubsealevelloadsvolume; - + return (sealevelloadsaverage+subsealevelloadsaverage)/totaloceanarea; } /*}}}*/ @@ -758,5 +752,5 @@ IssmDouble m1, m2, m3; bool rotation=false; - + /*early return?:*/ femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); @@ -786,8 +780,8 @@ ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - + ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - + ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -817,11 +811,10 @@ int* indices=xNew(nel); for(int i=0;i* vloadcopy=new Vector(nel); IssmDouble* loadcopy=xNew(nel); - + vloadcopy->SetValues(nel,indices,load,INS_VAL); vloadcopy->Assemble(); - if(subload){ Index: /issm/trunk-jpl/src/c/main/esmfbinders.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/esmfbinders.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/main/esmfbinders.cpp (revision 26468) @@ -115,5 +115,5 @@ Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); surface_input->GetInputAverage(&surface); - + *(issmoutputs+f*numberofelements+i) = surface; Index: /issm/trunk-jpl/src/c/main/issm.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/issm.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/main/issm.cpp (revision 26468) @@ -15,5 +15,5 @@ /*Initialize femmodel from arguments provided command line: */ FemModel *femmodel = new FemModel(argc,argv,comm_init); - + /*Need to know we are firing up from ISSM main, not a coupler driver like issm_slcp or issm_ocean:*/ femmodel->parameters->AddObject(new IntParam(IsSlcCouplingEnum,0)); Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26468) @@ -12,5 +12,5 @@ int melt_parameterization; femmodel->parameters->FindParam(&melt_parameterization,FrontalForcingsParamEnum); - + /*Calculate melting rate*/ switch(melt_parameterization){ Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp (revision 26468) @@ -27,6 +27,4 @@ } } - - void InputUpdateFromConstantx(FemModel* femmodel,IssmDouble constant, int name){ @@ -84,5 +82,4 @@ else _error_("InputUpdateFromConstantx error message: type not supported yet!"); - } void InputUpdateFromConstantx(Inputs* inputs,Elements* elements,bool constant, int name){ Index: /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp (revision 26468) @@ -49,5 +49,4 @@ descriptor=variables_descriptors[i]; - /*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */ if (strncmp(descriptor,"scaled_",7)==0){ @@ -166,5 +165,4 @@ _printf0_("\n"); } - if (femmodel->inputs->GetInputObjectEnum(MasstransportSpcthicknessEnum)==DatasetInputEnum) @@ -286,4 +284,3 @@ } - } /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/MmeToInputFromIdx/MmeToInputFromIdx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/MmeToInputFromIdx/MmeToInputFromIdx.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/MmeToInputFromIdx/MmeToInputFromIdx.cpp (revision 26468) @@ -7,5 +7,4 @@ #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif - #include "../../shared/shared.h" Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 26468) @@ -289,5 +289,4 @@ N_all[i] = N; - for(Object* & object : elements->objects){ Element* element=xDynamicCast(object); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 26468) @@ -717,5 +717,5 @@ xDelete(cost_functions_weights_N); xDelete(cost_functions_weights); - + /*}}}*/ } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26468) @@ -90,5 +90,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.issampling",TransientIssamplingEnum)); - + /*For stress balance only*/ parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum)); Index: /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 26468) @@ -24,5 +24,4 @@ int solutiontype; char* solutiontypestring = NULL; - /*recover my_rank:*/ @@ -71,5 +70,5 @@ char* root=NULL; char* modelname=NULL; - + femmodel->parameters->FindParam(&currEvalId,QmuCurrEvalIdEnum); femmodel->parameters->FindParam(&statistics,QmuStatisticsEnum); Index: /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp (revision 26468) @@ -20,5 +20,5 @@ IssmDouble* dmatfield = NULL; int* imatfield = NULL; - + //type of the returned field: int type; @@ -109,5 +109,5 @@ int range,lower_row,upper_row; int nfilesperdirectory; - + /*intermediary:*/ IssmDouble* doublemat=NULL; @@ -237,5 +237,5 @@ char* field=fields[f]; meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]); - + if(meanxtype[f]==1){ meanxsize[f]=1; @@ -246,5 +246,5 @@ timemean+=scalar/nsteps; } - + /*Figure out max and min of time means: */ if(i==(lower_row+1)){ @@ -334,5 +334,5 @@ IssmDouble* allmax=xNew(size); IssmDouble* allmin=xNew(size); - + ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); @@ -344,5 +344,5 @@ } /*}}}*/ } - + /*Now do the same for the time mean fields:*/ for (int f=0;f(size*nbins); - + ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); xDelete(histo); @@ -620,5 +619,5 @@ } _printf0_("Start aggregating time mean histogram:\n"); - + /*Now do the same for the time mean fields:*/ for (int f=0;f(size); IssmDouble* sumx2=xNew(size); - + ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); @@ -980,5 +979,5 @@ IssmDouble* sumx=xNew(size); IssmDouble* sumx2=xNew(size); - + ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); @@ -1021,5 +1020,5 @@ int* indices=NULL; int nindices; - + /*intermediary:*/ IssmDouble* doublemat=NULL; @@ -1143,5 +1142,5 @@ if(my_rank==0){ char fieldname[1000]; - + sprintf(fieldname,"%s%s",fields[f],"Samples"); results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allx,nsamples,1,j+1,0)); @@ -1152,7 +1151,7 @@ x=xs[counter]; allx=xNew(nsamples*nindices); - + MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm()); - + /*add to results:*/ if(my_rank==0){ @@ -1183,5 +1182,5 @@ FemModel* femmodel=new FemModel(); - + /*Some parameters that will allow us to use the OutputResultsx module:*/ parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false)); @@ -1327,8 +1326,8 @@ results = new Results(); parameters=new Parameters(); - + //solution type: parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum)); - + //root directory directory=xNew(strlen(argv[2])+1); @@ -1359,5 +1358,4 @@ ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm); - iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics"); for (int i=1;i<=numstatistics;i++){ @@ -1367,5 +1365,5 @@ int nsamples; _printf0_("Dealing with qmu statistical computation #" << i << "\n"); - + sprintf(string,"md.qmu.statistics.method(%i).name",i); iomodel->FindConstant(&name,string); @@ -1391,5 +1389,5 @@ iomodel->FetchData(&indices,&dummy,&nindices,string); parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices)); - + ComputeSampleSeries(parameters,results,color,statcomm); } @@ -1410,5 +1408,5 @@ /*output results:*/ OutputStatistics(parameters,results,color,statcomm); - + /*all meet here: */ ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n"); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 26468) @@ -1618,5 +1618,5 @@ //// MELT, PERCOLATION AND REFREEZE - + // initialize refreeze, runoff, flxDn and dW vectors [kg] IssmDouble* F = xNewZeroInit(n); @@ -1707,5 +1707,5 @@ for(int l=i;(l=dPHC-Dtol);l++) depthice=depthice+dz[l]; } - + // break loop if there is no meltwater and if depth is > mw_depth if (fabs(inM) < Wtol && i > X){ Index: /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalanceSicopolis.cpp (revision 26468) @@ -32,5 +32,5 @@ IssmDouble Pmax = 0.6; IssmDouble inv_twelve=1./12.; - + sconv=(rho_water/rho_ice); //rhow_rain/rhoi @@ -91,5 +91,5 @@ IssmDouble coeff5=3.8e-05; // polynomial IssmDouble coeff6=6.0e-07; // fit - + if(tstar>=temp_rain) frac_solid = 0.0; @@ -100,5 +100,5 @@ +tstar*(coeff3+tstar*(coeff4+tstar*(coeff5+tstar*coeff6)))); } - + snowfall=snowfall+precip_star*frac_solid*inv_twelve; } @@ -108,5 +108,5 @@ if(snowfall<0.0) snowfall=0.0; // correction of if(rainfall<0.0) rainfall=0.0; // negative values - + if(rainfall<=(Pmax*snowfall)){ if((rainfall+beta1*pdd)<=(Pmax*snowfall)) { @@ -126,5 +126,5 @@ runoff = smelt+rainfall-Pmax*snowfall; } - + saccu = precip; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp (revision 26468) @@ -76,5 +76,5 @@ Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters); femmodel->profiler->Stop(SOLVER); - + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys; @@ -98,5 +98,5 @@ } } - + /*Increase count: */ count++; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp (revision 26468) @@ -33,5 +33,4 @@ IssmPDouble t1,t2; /* Time measurement for bottleneck analysis */ - double tmp1,tmp2,tmp3; int tmpi; @@ -39,6 +38,4 @@ int noIt; - - int precond = 0; @@ -49,5 +46,4 @@ PetscBool flag,flg; #endif - char ksp_type[50]; @@ -71,5 +67,4 @@ #endif - if(precond){ _printf0_("Running WITH preconditioner\n"); @@ -78,8 +73,7 @@ } - /*Initialize output*/ Vector* out_uf=new Vector(uf0); - + /* Extract block matrices from the saddle point matrix */ /* [ A B ] = Kff @@ -97,5 +91,5 @@ MatGetSubMatrix(Kff,isp,isv,MAT_INITIAL_MATRIX,&BT); #endif - + /* Extract preconditioner matrix on the pressure space*/ #if PETSC_VERSION_GT(3,8,0) @@ -105,5 +99,4 @@ #endif - /* Get number of velocity / pressure nodes */ MatGetSize(B,&nu,&np); @@ -117,5 +110,4 @@ VecGetSubVector(out_uf->pvector->vector,isv,&uold); VecGetSubVector(out_uf->pvector->vector,isp,&p); - /* Set up intermediaries */ @@ -133,5 +125,5 @@ VecDuplicate(p,&wnew);VecSet(wnew,0.0); VecDuplicate(uold,&chi);VecSet(chi,0.0); - + /* Get global RHS (for each block sub-problem respectively)*/ VecGetSubVector(pf,isv,&f1); @@ -143,5 +135,5 @@ /* a(u0,v) = f1(v)-b(p0,v) i.e. Au0 = F1-Bp0 */ /* u0 = u_DIR on \Gamma_DIR */ - + /* Create KSP context */ KSPCreate(IssmComm::GetComm(),&kspu); @@ -174,13 +166,11 @@ _error_("Suggested KSP method not implemented yet!\n"); } - + KSPSetInitialGuessNonzero(kspu,PETSC_TRUE); - + /*Strong rel. residual tolerance needed when solving for the velocity update. * This is because ISSM uses the dimensional equations, so the initial guess chi = 0 * results in a very high initial residual (~ 1e19)*/ KSPSetTolerances(kspu,ELLTOL,PETSC_DEFAULT,PETSC_DEFAULT,maxiter); //maxiter - - KSPGetPC(kspu,&pcu); @@ -210,5 +200,4 @@ KSPSetUp(kspu); - /* Create RHS */ /* RHS = F1-B * pold */ @@ -227,7 +216,5 @@ VecNorm(resu,NORM_2,&tmp5); - - } - + } /* Go solve Au0 = F1-Bp0*/ @@ -235,9 +222,7 @@ KSPGetIterationNumber(kspu,&noIt); - - if (VerboseConvergence()) { - + KSPGetIterationNumber(kspu,&tmpi); VecScale(rhsu,-1.);MatMultAdd(A,uold,rhsu,tmpu);VecScale(rhsu,-1.); @@ -249,17 +234,13 @@ } - - /* Set up u_new */ VecDuplicate(uold,&unew);VecCopy(uold,unew); VecAssemblyBegin(unew);VecAssemblyEnd(unew); - - /* ------------------------------------------------------------- */ /*Get initial residual*/ /*(1/mu(x) * g0, q) = b(q,u0) - (f2,q) i.e. IP * g0 = BT * u0 - F2*/ - + /* Create KSP context */ KSPCreate(IssmComm::GetComm(),&kspip); @@ -269,15 +250,13 @@ KSPSetOperators(kspip,IP,IP,DIFFERENT_NONZERO_PATTERN); #endif - + /* Create RHS */ /* RHS = BT * uold - F2 */ VecScale(uold,-1.);MatMultAdd(BT,uold,f2,rhsp);VecScale(uold,-1.); - /* Set KSP & PC options */ KSPSetType(kspip,KSPGMRES); KSPSetInitialGuessNonzero(kspip,PETSC_TRUE); - - + KSPGetPC(kspip,&pcp); PCSetType(pcp,PCJACOBI); @@ -294,14 +273,10 @@ } - - /* Go solve */ KSPSolve(kspip,rhsp,gold); - - if (VerboseConvergence()) { - + KSPGetResidualNorm(kspip,&tmp1); VecScale(rhsp,-1.);MatMultAdd(IP,gold,rhsp,tmpp);VecScale(rhsp,-1.); @@ -311,5 +286,4 @@ _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: "< 1e-5) @@ -326,5 +299,4 @@ VecDuplicate(gold,&gnew);VecCopy(gold,gnew); VecAssemblyBegin(gnew);VecAssemblyEnd(gnew); - /* ------------------------------------------------------------ */ @@ -334,5 +306,4 @@ VecDuplicate(gold,&wold);VecCopy(gold,wold); VecAssemblyBegin(wold);VecAssemblyEnd(wold); - /* Count number of iterations */ @@ -351,11 +322,10 @@ VecScale(wold,1.);MatMult(B,wold,rhsu);VecScale(wold,1.); VecSet(chi,0.0); - - + if(VerboseConvergence()) { VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.); VecNorm(tmpu,NORM_2,&tmp4); - + KSPInitialResidual(kspu,chi,tmpu,tmpu2,resu,rhsu); VecNorm(resu,NORM_2,&tmp5); @@ -363,26 +333,19 @@ } - KSPSolve(kspu,rhsu,chi); - - - - + if (VerboseConvergence()) { VecNorm(chi,NORM_2,&tmp1); KSPGetResidualNorm(kspu,&tmp2); - + VecScale(rhsu,-1.);MatMultAdd(A,chi,rhsu,tmpu);VecScale(rhsu,-1.); VecNorm(tmpu,NORM_2,&tmp3); - - + KSPGetIterationNumber(kspu,&tmpi); _printf0_("||chi_nk||_euc: "<< tmp1 << "\n||A*chi_0k - rhsu||_euc: "< 10./TOL) @@ -461,8 +416,6 @@ break; } - /* ---------------------------------------------------------- */ - /*Directional update*/ @@ -475,5 +428,5 @@ /* Assign new to old iterates */ VecCopy(wnew,wold);VecCopy(gnew,gold);VecCopy(unew,uold); - + count++; } @@ -487,10 +440,9 @@ *puf=out_uf; - /* Cleanup */ KSPDestroy(&kspu);KSPDestroy(&kspip); MatDestroy(&A);MatDestroy(&B);MatDestroy(&BT);MatDestroy(&IP); - + VecDestroy(&p);VecDestroy(&uold);VecDestroy(&unew);VecDestroy(&rhsu);VecDestroy(&rhsp); VecDestroy(&gold);VecDestroy(&gnew);VecDestroy(&wold);VecDestroy(&wnew);VecDestroy(&chi); @@ -522,5 +474,4 @@ double rnorm1, rnorm2; - if(VerboseModule()) _printf0_(" checking convergence\n"); @@ -548,5 +499,5 @@ MatGetSubMatrix(Kff->pmatrix->matrix,isp,isv,MAT_INITIAL_MATRIX,&BT); #endif - + /*no. of free nodes in velocity/pressure space*/ MatGetSize(B,&n_u,&n_p); @@ -560,5 +511,4 @@ VecGetSubVector(uf->pvector->vector,isv,&u); VecGetSubVector(uf->pvector->vector,isp,&p); - /*Extract values of the RHS corresponding to the first/second block*/ @@ -573,13 +523,11 @@ VecDuplicate(p,&res2);VecSet(res2,1.0); - /*Display solver caracteristics*/ if (VerboseConvergence()){ - + /*Calculate res1 = A*u + B*p - f1*/ VecScale(f1,-1.);MatMultAdd(A,u,f1,tmp);MatMultAdd(B,p,tmp,res1);VecScale(f1,-1.); /*Calculate res2 = B^T * u - f2*/ VecScale(f2,-1.);MatMultAdd(BT,u,f2,res2);VecScale(f2,-1.); - /*compute norm(res1), norm(res2), norm(F) and residue*/ @@ -597,9 +545,8 @@ VecRestoreSubVector(uf->pvector->vector,isv,&u); VecRestoreSubVector(uf->pvector->vector,isp,&p); - + /*Extract values corresponding to velocity/pressure on the old solution*/ VecGetSubVector(old_uf->pvector->vector,isv,&uold); VecGetSubVector(old_uf->pvector->vector,isp,&pold); - /*Force equilibrium (Mandatory)*/ @@ -609,5 +556,5 @@ /*Calculate res2 = B^T * uold - f2*/ VecScale(f2,-1.);MatMultAdd(BT,uold,f2,res2);VecScale(f2,-1.); - + /*compute norm(res1), norm(res2), norm(F) and residue*/ VecNorm(res1,NORM_2,&rnorm1);VecNorm(res2,NORM_2,&rnorm2); @@ -626,6 +573,4 @@ VecRestoreSubVector(old_uf->pvector->vector,isv,&uold); VecRestoreSubVector(old_uf->pvector->vector,isp,&pold); - - //print @@ -721,5 +666,5 @@ int count=0; bool converged=false; - + /*Start non-linear iteration using input velocity: */ GetSolutionFromInputsx(&ug,femmodel); @@ -729,5 +674,5 @@ InputUpdateFromConstantx(femmodel,converged,ConvergedEnum); InputUpdateFromSolutionx(femmodel,ug); - + for(;;){ @@ -747,7 +692,7 @@ PetscOptionsGetInt(NULL,PETSC_NULL,"-schur_pc",&precond,NULL); #endif - + StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); - + /* Construct Schur preconditioner matrix or mass matrix depending on input */ if(precond) @@ -760,5 +705,5 @@ } }else{ - + for(Object* & object : femmodel->elements->objects){ Element* element2=xDynamicCast(object); @@ -769,7 +714,6 @@ } - delete analysis; - + /*Obtain index sets for velocity and pressure components */ IS isv = NULL; @@ -785,9 +729,7 @@ Kff->Assemble(); - /*Solve*/ femmodel->profiler->Start(SOLVER); _assert_(Kff->type==PetscMatType); - SchurCGSolver(&uf, Index: /issm/trunk-jpl/src/c/toolkits/codipack/ampi_interface.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/codipack/ampi_interface.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/toolkits/codipack/ampi_interface.cpp (revision 26468) @@ -12,3 +12,2 @@ #error "Cannot compile without MeDiPack and AdjointMpi" #endif - Index: /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp (revision 26467) +++ /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp (revision 26468) @@ -318,5 +318,5 @@ xDelete(valueX); xDelete(indexX); - + delete data; }