Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/MasstransportAnalysis.cpp (revision 22142) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/MasstransportAnalysis.cpp (revision 22143) @@ -157,4 +157,15 @@ iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum); + + /*Initialize cumdeltalthickness input: unfortunately, we don't have femmodel, so we need to iterate on + *elements, otherwise we would have used InputUpdateFromConstantx: */ + counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum); + counter++; + } + } if(!issmb){ @@ -650,5 +661,6 @@ basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); IssmDouble* newthickness = xNew(numnodes); - IssmDouble* deltathickness = xNew(numnodes); + IssmDouble* cumdeltathickness = xNew(numnodes); + IssmDouble* deltathickness = xNew(numnodes); IssmDouble* newbase = xNew(numnodes); IssmDouble* newsurface = xNew(numnodes); @@ -675,7 +687,11 @@ basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum); basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum); - - /*What is the delta thickness forcing the sea-level rise core: */ - for(i=0;iGetInputListOnNodes(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum); + + /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/ + for(i=0;iAddBasalInput(ThicknessEnum,newthickness,P1Enum); - element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum); element->AddBasalInput(SurfaceEnum,newsurface,P1Enum); element->AddBasalInput(BaseEnum,newbase,P1Enum); + element->AddBasalInput(SealevelriseCumDeltathicknessEnum,cumdeltathickness,P1Enum); + element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum); /*Free ressources:*/ Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp (revision 22142) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp (revision 22143) @@ -39,5 +39,5 @@ iomodel->FetchDataToInput(elements,"md.slr.spcthickness",SealevelriseSpcthicknessEnum); iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0); - + /*Initialize cumdeltalthickness input: unfortunately, we don't have femmodel, so we need to iterate on *elements, otherwise we would have used InputUpdateFromConstantx: */ @@ -88,5 +88,4 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.slr.run_frequency",SealevelriseRunFrequencyEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum)); @@ -97,5 +96,4 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum)); - parameters->SetParam(1,SealevelriseRunCountEnum); iomodel->FetchData(&elastic,"md.slr.elastic"); Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h (revision 22142) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h (revision 22143) @@ -66,5 +66,4 @@ void TransferSealevel(FemModel* femmodel,int forcingenum); void EarthMassTransport(FemModel* femmodel); -bool SlrAccumulateDeltaThickness(FemModel* femmodel); //solution configuration Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp (revision 22142) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp (revision 22143) @@ -19,4 +19,5 @@ Vector *U_north = NULL; Vector *U_east = NULL; + Vector* cumdeltathickness=NULL; bool save_results,isslr,iscoupler; int configuration_type; @@ -35,5 +36,5 @@ IssmDouble *zz = NULL; IssmDouble dt,steric_rate; - int frequency; + int frequency,count; int loop; int horiz; @@ -45,4 +46,7 @@ femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); + femmodel->parameters->FindParam(&frequency,SealevelriseRunFrequencyEnum); + femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum); + /*first, recover lat,long and radius vectors from vertices: */ @@ -77,5 +81,5 @@ if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); - /*figure out whether deltathickness was computed on the earth (mass trasnport module should + /*figure out whether deltathickness was computed on the earth (mass transport module should * have taken care of it for ice caps: */ if(iscoupler){ @@ -88,9 +92,18 @@ } - /*transfer deltathickness forcing from ice caps to earth model: */ - if(iscoupler) TransferForcing(femmodel,SealevelriseDeltathicknessEnum); + /*transfer cum deltathickness forcing from ice caps to earth model: */ + if(iscoupler & (count==frequency)) TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum); /*call sea-level rise sub cores:*/ - if(isslr && !SlrAccumulateDeltaThickness(femmodel)){ + if(isslr & (count==frequency)){ + + if(VerboseSolution()) _printf0_(" computing SLR core\n"); + + if(iscoupler){ + /*we have accumulated thicknesses, dump them in deltathcikness: */ + GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum); + delete cumdeltathickness; + } femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); @@ -154,5 +167,27 @@ } - + + /*transfer sea-level back to ice caps, reset cum thicknesses to 0 and reset counter: */ + if(iscoupler){ + if (count==frequency) { + //transfer sea level back to ice caps: + TransferSealevel(femmodel,SealevelEnum); + + //reset cumdeltathickness to 0: + GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); + cumdeltathickness->Set(0); + InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); + delete cumdeltathickness; + + //reset counter to 1: + femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter. + } + else{ + //increment counter: + count++; + femmodel->parameters->SetParam(count,SealevelriseRunCountEnum); + } + } + /*Free ressources: */ delete longitude; @@ -163,6 +198,4 @@ delete zz; - /*transfer sea-level back to ice caps: */ - if(iscoupler)TransferSealevel(femmodel,SealevelEnum); } /*}}}*/ @@ -260,5 +293,5 @@ /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular * ice cap: */ - forcingglobal->SetValues(M,index,forcingfromcap,INS_VAL); + forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL); xDelete(index); } @@ -406,4 +439,5 @@ Vector *newthickness = NULL; Vector *deltathickness = NULL; + Vector *cumdeltathickness = NULL; int nv; @@ -430,61 +464,14 @@ InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum); + /*add to cumulated delta thickness: */ + GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); + cumdeltathickness->AXPY(deltathickness,1); + InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); + /*free ressources:*/ delete oldthickness; delete newthickness; delete deltathickness; + delete cumdeltathickness; } /*}}}*/ -bool SlrAccumulateDeltaThickness(FemModel* femmodel){ /*{{{*/ - - Vector* deltathickness=NULL; - Vector* cumdeltathickness=NULL; - int frequency=1; - int count; - int solution_type; - - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - if(solution_type==SealevelriseSolutionEnum)return false ; //we don't accumulate when we run this solution type! - - femmodel->parameters->FindParam(&frequency,SealevelriseRunFrequencyEnum); - femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum); - - if (count==frequency){ - femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter. - - /*increment delta thickness load: */ - GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); - GetVectorFromInputsx(&deltathickness,femmodel,SealevelriseDeltathicknessEnum,VertexSIdEnum); - - cumdeltathickness->AXPY(deltathickness,1); - - /*update deltathickness with the final cumulated thickness: */ - InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum); - - /*reset cumulative thickness to 0: */ - cumdeltathickness->Set(0); - InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); - - delete deltathickness; - delete cumdeltathickness; - - return false; - } - else{ - count++; - femmodel->parameters->SetParam(count,SealevelriseRunCountEnum); //increment counter - - /*increment cumulative delta thickness load: */ - GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); - GetVectorFromInputsx(&deltathickness,femmodel,SealevelriseDeltathicknessEnum,VertexSIdEnum); - - cumdeltathickness->AXPY(deltathickness,1); - InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); - - delete deltathickness; - delete cumdeltathickness; - - return true; - } - -} /*}}}*/ Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22142) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22143) @@ -66,4 +66,6 @@ parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.run_frequency",SealevelriseRunFrequencyEnum)); + parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); {/*This is specific to ice...*/ parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum)); Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m (revision 22142) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m (revision 22143) @@ -58,4 +58,13 @@ end end + + %check that run_frequency is the same everywhere: + for i=1:length(slm.icecaps), + if slm.icecaps{i}.slr.run_frequency~=slm.earth.slr.run_frequency, + error(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name)); + end + end + +