Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 20135) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 20136) @@ -51,4 +51,10 @@ int numoutputs; char** requestedoutputs = NULL; + + /*transition vectors: */ + IssmDouble **transitions = NULL; + int *transitions_M = NULL; + int *transitions_N = NULL; + int ntransitions; /*some constant parameters: */ @@ -126,4 +132,18 @@ } + /*Transitions: */ + iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,SealevelriseTransitionsEnum); + if(transitions){ + parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N)); + + for(int i=0;i(transition); + } + xDelete(transitions); + xDelete(transitions_M); + xDelete(transitions_N); + } + /*Requested outputs*/ iomodel->FetchData(&requestedoutputs,&numoutputs,SealevelriseRequestedOutputsEnum); Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 20135) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 20136) @@ -50,6 +50,6 @@ void damage_core(FemModel* femmodel); void sealevelrise_core(FemModel* femmodel); -Vector * sealevelrise_core_eustatic(FemModel* femmodel); -Vector * sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic); +Vector* sealevelrise_core_eustatic(FemModel* femmodel); +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic); IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel); @@ -62,4 +62,6 @@ void ResetBoundaryConditions(FemModel* femmodel, int analysis_type); void PrintBanner(void); +void TransferForcing(FemModel* femmodel,int forcingenum); +void TransferSealevel(FemModel* femmodel,int forcingenum); //solution configuration Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 20135) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 20136) @@ -10,42 +10,280 @@ #include "../solutionsequences/solutionsequences.h" -void sealevelrise_core(FemModel* femmodel){ +void sealevelrise_core(FemModel* femmodel){ /*{{{*/ Vector *Sg = NULL; Vector *Sg_eustatic = NULL; - bool save_results; + bool save_results,isslr,iscoupler; int configuration_type; + int solution_type; int numoutputs = 0; char **requested_outputs = NULL; - - if(VerboseSolution()) _printf0_(" computing sea level rise\n"); - + /*Recover some parameters: */ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - femmodel->parameters->FindParam(&numoutputs,SealevelriseNumRequestedOutputsEnum); - if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum); + femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); + femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); + + /*several cases here, depending on value of iscoupler and isslr: + ( !iscoupler & !isslr) we are not interested in being here :) + ( !iscoupler & isslr) we are running in uncoupled mode + ( iscoupler & isslr) we are running in coupled mode, and better be earth + ( iscoupler & !isslr) we are running in coupled mode, and better be an ice cap + */ + + /*early return: */ + if( !iscoupler & !isslr) return; //we are not interested in being here :) + + /*In what follows we assume we are all running slr, either in coupled, or uncoupled mode:*/ + if(VerboseSolution()) _printf0_(" computing sea level rise\n"); /*set configuration: */ - femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); - - /*call two sub cores:*/ - Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS. - - Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //sea-level rise dependent terms (2nd and 5th terms on the RHS) - - /*get results out:*/ - InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum); - - if(save_results){ - if(VerboseSolution()) _printf0_(" saving results\n"); - femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); - } - - /*Free ressources:*/ - if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} - delete Sg; - delete Sg_eustatic; - - -} + if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); + + /*transfer deltathickness forcing from ice caps to earth model: */ + if(iscoupler) TransferForcing(femmodel,SealevelriseDeltathicknessEnum); + + /*call sea-level rise sub cores:*/ + if(isslr){ + Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS. + + Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //sea-level rise dependent terms (2nd and 5th terms on the RHS) + + /*get results into elements:*/ + InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + } + + if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + delete Sg; + delete Sg_eustatic; + if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} + } + + /*transfer sea-level back to ice caps: */ + if(iscoupler)TransferSealevel(femmodel,SealevelriseSEnum); +} +/*}}}*/ +void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/ + + /*forcing being transferred from models to earth: */ + IssmDouble** forcings=NULL; + IssmDouble* forcing=NULL; + Vector* forcingglobal=NULL; + int* nvs=NULL; + + /*transition vectors:*/ + IssmDouble** transitions=NULL; + int ntransitions; + int* transitions_m=NULL; + int* transitions_n=NULL; + int nv; + + /*communicators:*/ + ISSM_MPI_Comm tocomm; + ISSM_MPI_Comm* fromcomms=NULL; + ISSM_MPI_Status status; + int my_rank; + int modelid,earthid; + int nummodels; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&modelid,ModelIdEnum); + femmodel->parameters->FindParam(&earthid,EarthIdEnum); + femmodel->parameters->FindParam(&nummodels,NumModelsEnum); + my_rank=IssmComm::GetRank(); + + /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */ + if(modelid==earthid)femmodel->parameters->FindParam(&fromcomms,&nv,IcecapToEarthCommEnum); + else femmodel->parameters->FindParam(&tocomm, IcecapToEarthCommEnum); + + /*For each icecap, retrieve the forcing vector that will be sent to the earth model: */ + if(modelid!=earthid){ + nv=femmodel->vertices->NumberOfVertices(); + GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum); + } + + /*Send the forcing to the earth model:{{{*/ + if(my_rank==0){ + if(modelid==earthid){ + forcings=xNew(nummodels-1); + nvs=xNew(nummodels-1); + for(int i=0;i(nvs[i]); + MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status); + } + + } + else{ + ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm); + ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, tocomm); + } + } + /*}}}*/ + + /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/ + if(modelid==earthid){ + + /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions. + *First, build the global delta thickness vector in the earth model: */ + nv=femmodel->vertices->NumberOfVertices(); + forcingglobal= new Vector(nv); + + /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/ + femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum); + + if(ntransitions!=earthid)_error_("TransferForcing error message: number of transition vectors is not equal to the number of icecaps!"); + + /*Go through all the delta thicknesses coming from each ice cap: */ + if(my_rank==0){ + for(int i=0;i(M); for(int i=0;i(transition[i])-1; //matlab indexing! + + + /*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); + xDelete(index); + } + } + + /*Assemble vector:*/ + forcingglobal->Assemble(); + + /*Plug into elements:*/ + InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum); + } + /*}}}*/ + + /*Free ressources:{{{*/ + if(forcings){ + for(int i=0;i(temp); + } + xDelete(forcings); + } + if(forcing)xDelete(forcing); + if(forcingglobal)delete forcingglobal; + if(transitions){ + for(int i=0;i(temp); + } + xDelete(transitions); + xDelete(transitions_m); + xDelete(transitions_n); + } + if(nvs)xDelete(nvs); + /*}}}*/ + +} /*}}}*/ +void TransferSealevel(FemModel* femmodel,int forcingenum){ /*{{{*/ + + /*forcing being transferred from earth to ice caps: */ + IssmDouble* forcing=NULL; + IssmDouble* forcingglobal=NULL; + + /*transition vectors:*/ + IssmDouble** transitions=NULL; + int ntransitions; + int* transitions_m=NULL; + int* transitions_n=NULL; + int nv; + + /*communicators:*/ + ISSM_MPI_Comm fromcomm; + ISSM_MPI_Comm* tocomms=NULL; + ISSM_MPI_Status status; + int my_rank; + int modelid,earthid; + int nummodels; + int numcoms; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&modelid,ModelIdEnum); + femmodel->parameters->FindParam(&earthid,EarthIdEnum); + femmodel->parameters->FindParam(&nummodels,NumModelsEnum); + my_rank=IssmComm::GetRank(); + + /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/ + if(modelid==earthid)femmodel->parameters->FindParam(&tocomms,&numcoms,IcecapToEarthCommEnum); + else femmodel->parameters->FindParam(&fromcomm, IcecapToEarthCommEnum); + + + /*Retrieve sea-level on earth model: */ + if(modelid==earthid){ + nv=femmodel->vertices->NumberOfVertices(); + GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum); + } + + /*Send the forcing to the ice caps:{{{*/ + if(my_rank==0){ + + if(modelid==earthid){ + + /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */ + femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum); + + if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!"); + + for(int i=0;i(nv); + IssmDouble* transition=transitions[i]; + for(int j=0;j(transition[j])-1]; + } + ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, i, tocomms[i]); + ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, i, tocomms[i]); + } + } + else{ + ISSM_MPI_Recv(&nv, 1, ISSM_MPI_INT, 0, modelid, fromcomm, &status); + forcing=xNew(nv); + ISSM_MPI_Recv(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, fromcomm, &status); + } + } + /*}}}*/ + + /*On each ice cap, spread the forcing across cpus, and update the elements dataset accordingly: {{{*/ + if(modelid!=earthid){ + + ISSM_MPI_Bcast(&nv,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + if(my_rank!=0)forcing=xNew(nv); + ISSM_MPI_Bcast(forcing,nv,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Plug into elements:*/ + InputUpdateFromVectorx(femmodel,forcing,forcingenum,VertexSIdEnum); + } + /*}}}*/ + + /*Free ressources:{{{*/ + if(forcingglobal)xDelete(forcingglobal); + if(forcing)xDelete(forcing); + if(transitions){ + for(int i=0;i(temp); + } + xDelete(transitions); + xDelete(transitions_m); + xDelete(transitions_n); + } + /*}}}*/ + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 20135) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 20136) @@ -21,5 +21,5 @@ /*parameters: */ IssmDouble starttime,finaltime,dt,yts; - bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalving; + bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isslr,iscoupler,islevelset,isdamageevolution,ishydrology,iscalving; bool save_results,dakota_analysis; bool time_adapt; @@ -51,4 +51,6 @@ femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum); femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum); + femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); + femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum); @@ -141,4 +143,7 @@ } + /*Sea level rise: */ + if(isslr | iscoupler) sealevelrise_core(femmodel); + /*unload results*/ if(VerboseSolution()) _printf0_(" computing requested outputs\n");