| [19984] | 1 | #include "./SealevelriseAnalysis.h"
|
|---|
| 2 | #include "../toolkits/toolkits.h"
|
|---|
| 3 | #include "../classes/classes.h"
|
|---|
| 4 | #include "../shared/shared.h"
|
|---|
| 5 | #include "../modules/modules.h"
|
|---|
| 6 |
|
|---|
| 7 | /*Model processing*/
|
|---|
| 8 | void SealevelriseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
|
|---|
| 9 | /*No constraints*/
|
|---|
| 10 | }/*}}}*/
|
|---|
| 11 | void SealevelriseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
|
|---|
| 12 | /*No loads*/
|
|---|
| 13 | }/*}}}*/
|
|---|
| 14 | void SealevelriseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
|
|---|
| 15 | ::CreateNodes(nodes,iomodel,SealevelriseAnalysisEnum,P1Enum);
|
|---|
| 16 | }/*}}}*/
|
|---|
| 17 | int SealevelriseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
|
|---|
| 18 | return 1;
|
|---|
| 19 | }/*}}}*/
|
|---|
| 20 | void SealevelriseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
|
|---|
| 21 |
|
|---|
| 22 | /*Update elements: */
|
|---|
| 23 | int counter=0;
|
|---|
| 24 | for(int i=0;i<iomodel->numberofelements;i++){
|
|---|
| 25 | if(iomodel->my_elements[i]){
|
|---|
| 26 | Element* element=(Element*)elements->GetObjectByOffset(counter);
|
|---|
| 27 | element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
|
|---|
| 28 | counter++;
|
|---|
| 29 | }
|
|---|
| 30 | }
|
|---|
| 31 |
|
|---|
| 32 | /*Create inputs: */
|
|---|
| [20690] | 33 | iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
|
|---|
| 34 | iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
|
|---|
| 35 | iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum);
|
|---|
| [20704] | 36 | iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum);
|
|---|
| [20690] | 37 | iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
|
|---|
| [21752] | 38 | iomodel->FetchDataToInput(elements,"md.slr.steric_rate",SealevelriseStericRateEnum);
|
|---|
| [19984] | 39 |
|
|---|
| 40 | }/*}}}*/
|
|---|
| 41 | void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
|
|---|
| 42 |
|
|---|
| [19986] | 43 | int nl;
|
|---|
| [19984] | 44 | IssmDouble* love_h=NULL;
|
|---|
| 45 | IssmDouble* love_k=NULL;
|
|---|
| [20999] | 46 | IssmDouble* love_l=NULL;
|
|---|
| [19986] | 47 |
|
|---|
| [20029] | 48 | bool elastic=false;
|
|---|
| [20028] | 49 | IssmDouble* G_elastic = NULL;
|
|---|
| [20033] | 50 | IssmDouble* G_elastic_local = NULL;
|
|---|
| [20999] | 51 | IssmDouble* U_elastic = NULL;
|
|---|
| 52 | IssmDouble* U_elastic_local = NULL;
|
|---|
| 53 | IssmDouble* H_elastic = NULL;
|
|---|
| 54 | IssmDouble* H_elastic_local = NULL;
|
|---|
| [20033] | 55 | int M,m,lower_row,upper_row;
|
|---|
| 56 | IssmDouble degacc=.01;
|
|---|
| [19984] | 57 |
|
|---|
| [20036] | 58 | int numoutputs;
|
|---|
| 59 | char** requestedoutputs = NULL;
|
|---|
| 60 |
|
|---|
| [20136] | 61 | /*transition vectors: */
|
|---|
| 62 | IssmDouble **transitions = NULL;
|
|---|
| 63 | int *transitions_M = NULL;
|
|---|
| 64 | int *transitions_N = NULL;
|
|---|
| 65 | int ntransitions;
|
|---|
| 66 |
|
|---|
| [19986] | 67 | /*some constant parameters: */
|
|---|
| [20690] | 68 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.reltol",SealevelriseReltolEnum));
|
|---|
| 69 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.abstol",SealevelriseAbstolEnum));
|
|---|
| 70 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.maxiter",SealevelriseMaxiterEnum));
|
|---|
| 71 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
|
|---|
| 72 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
|
|---|
| [21331] | 73 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
|
|---|
| 74 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_h",SealevelriseTidalLoveHEnum));
|
|---|
| 75 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.tide_love_k",SealevelriseTidalLoveKEnum));
|
|---|
| [21344] | 76 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.fluid_love",SealevelriseFluidLoveEnum));
|
|---|
| 77 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.equatorial_moi",SealevelriseEquatorialMoiEnum));
|
|---|
| 78 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.polar_moi",SealevelrisePolarMoiEnum));
|
|---|
| [21331] | 79 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
|
|---|
| [21295] | 80 | parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
|
|---|
| [19986] | 81 |
|
|---|
| [20690] | 82 | iomodel->FetchData(&elastic,"md.slr.elastic");
|
|---|
| [20029] | 83 | if(elastic){
|
|---|
| [19984] | 84 |
|
|---|
| [20029] | 85 | /*love numbers: */
|
|---|
| [20690] | 86 | iomodel->FetchData(&love_h,&nl,NULL,"md.slr.love_h");
|
|---|
| 87 | iomodel->FetchData(&love_k,&nl,NULL,"md.slr.love_k");
|
|---|
| [20999] | 88 | iomodel->FetchData(&love_l,&nl,NULL,"md.slr.love_l");
|
|---|
| [20022] | 89 |
|
|---|
| [20029] | 90 | /*compute elastic green function for a range of angles*/
|
|---|
| [20690] | 91 | iomodel->FetchData(°acc,"md.slr.degacc");
|
|---|
| [20034] | 92 | M=reCast<int,IssmDouble>(180./degacc+1.);
|
|---|
| [21742] | 93 |
|
|---|
| 94 | // AD performance is sensitive to calls to ensurecontiguous.
|
|---|
| 95 | // // Providing "t" will cause ensurecontiguous to be called.
|
|---|
| 96 | #ifdef _HAVE_AD_
|
|---|
| [21741] | 97 | G_elastic=xNew<IssmDouble>(M,"t");
|
|---|
| 98 | U_elastic=xNew<IssmDouble>(M,"t");
|
|---|
| 99 | H_elastic=xNew<IssmDouble>(M,"t");
|
|---|
| [21742] | 100 | #else
|
|---|
| 101 | G_elastic=xNew<IssmDouble>(M);
|
|---|
| 102 | U_elastic=xNew<IssmDouble>(M);
|
|---|
| 103 | H_elastic=xNew<IssmDouble>(M);
|
|---|
| 104 | #endif
|
|---|
| 105 |
|
|---|
| [20044] | 106 |
|
|---|
| 107 | /*compute combined legendre + love number (elastic green function:*/
|
|---|
| [20033] | 108 | m=DetermineLocalSize(M,IssmComm::GetComm());
|
|---|
| 109 | GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
|
|---|
| [21742] | 110 | // AD performance is sensitive to calls to ensurecontiguous.
|
|---|
| 111 | // // Providing "t" will cause ensurecontiguous to be called.
|
|---|
| 112 | #ifdef _HAVE_AD_
|
|---|
| [21741] | 113 | G_elastic_local=xNew<IssmDouble>(m,"t");
|
|---|
| 114 | U_elastic_local=xNew<IssmDouble>(m,"t");
|
|---|
| 115 | H_elastic_local=xNew<IssmDouble>(m,"t");
|
|---|
| [21742] | 116 | #else
|
|---|
| 117 | G_elastic_local=xNew<IssmDouble>(m);
|
|---|
| 118 | U_elastic_local=xNew<IssmDouble>(m);
|
|---|
| 119 | H_elastic_local=xNew<IssmDouble>(m);
|
|---|
| 120 | #endif
|
|---|
| [20033] | 121 |
|
|---|
| 122 | for(int i=lower_row;i<upper_row;i++){
|
|---|
| [20029] | 123 | IssmDouble alpha,x;
|
|---|
| 124 | alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
|
|---|
| 125 |
|
|---|
| [20033] | 126 | G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
|
|---|
| [20999] | 127 | U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0);
|
|---|
| 128 | H_elastic_local[i-lower_row]= 0;
|
|---|
| [21908] | 129 | IssmDouble Pn = 0.;
|
|---|
| 130 | IssmDouble Pn1 = 0.;
|
|---|
| 131 | IssmDouble Pn2 = 0.;
|
|---|
| 132 | IssmDouble Pn_p = 0.;
|
|---|
| 133 | IssmDouble Pn_p1 = 0.;
|
|---|
| 134 | IssmDouble Pn_p2 = 0.;
|
|---|
| 135 |
|
|---|
| [20033] | 136 | for (int n=0;n<nl;n++) {
|
|---|
| [20999] | 137 | IssmDouble deltalove_G;
|
|---|
| 138 | IssmDouble deltalove_U;
|
|---|
| [20033] | 139 |
|
|---|
| [20999] | 140 | deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
|
|---|
| 141 | deltalove_U = (love_h[n]-love_h[nl-1]);
|
|---|
| 142 |
|
|---|
| 143 | /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
|
|---|
| 144 | if(n==0){
|
|---|
| 145 | Pn=1;
|
|---|
| 146 | Pn_p=0;
|
|---|
| 147 | }
|
|---|
| 148 | else if(n==1){
|
|---|
| 149 | Pn = cos(alpha);
|
|---|
| 150 | Pn_p = 1;
|
|---|
| 151 | }
|
|---|
| 152 | else{
|
|---|
| 153 | Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
|
|---|
| 154 | Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
|
|---|
| 155 | }
|
|---|
| [20033] | 156 | Pn2=Pn1; Pn1=Pn;
|
|---|
| [20999] | 157 | Pn_p2=Pn_p1; Pn_p1=Pn_p;
|
|---|
| [20033] | 158 |
|
|---|
| [20999] | 159 | G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
|
|---|
| 160 | U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
|
|---|
| 161 | H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
|
|---|
| [20033] | 162 | }
|
|---|
| 163 | }
|
|---|
| 164 |
|
|---|
| [20999] | 165 | /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
|
|---|
| [20033] | 166 | int* recvcounts=xNew<int>(IssmComm::GetSize());
|
|---|
| 167 | int* displs=xNew<int>(IssmComm::GetSize());
|
|---|
| 168 |
|
|---|
| 169 | //recvcounts:
|
|---|
| 170 | ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
|
|---|
| 171 |
|
|---|
| 172 | /*displs: */
|
|---|
| 173 | ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
|
|---|
| 174 |
|
|---|
| 175 | /*All gather:*/
|
|---|
| 176 | ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
|
|---|
| [20999] | 177 | ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
|
|---|
| 178 | ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
|
|---|
| [20033] | 179 | /*free ressources: */
|
|---|
| 180 | xDelete<int>(recvcounts);
|
|---|
| 181 | xDelete<int>(displs);
|
|---|
| 182 |
|
|---|
| 183 | /*}}}*/
|
|---|
| 184 |
|
|---|
| [20029] | 185 | /*Avoid singularity at 0: */
|
|---|
| 186 | G_elastic[0]=G_elastic[1];
|
|---|
| 187 | parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
|
|---|
| [20999] | 188 | U_elastic[0]=U_elastic[1];
|
|---|
| 189 | parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
|
|---|
| 190 | H_elastic[0]=H_elastic[1];
|
|---|
| 191 | parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
|
|---|
| [20028] | 192 |
|
|---|
| [20029] | 193 | /*free ressources: */
|
|---|
| 194 | xDelete<IssmDouble>(love_h);
|
|---|
| 195 | xDelete<IssmDouble>(love_k);
|
|---|
| [20999] | 196 | xDelete<IssmDouble>(love_l);
|
|---|
| [20029] | 197 | xDelete<IssmDouble>(G_elastic);
|
|---|
| [20033] | 198 | xDelete<IssmDouble>(G_elastic_local);
|
|---|
| [20999] | 199 | xDelete<IssmDouble>(U_elastic);
|
|---|
| 200 | xDelete<IssmDouble>(U_elastic_local);
|
|---|
| 201 | xDelete<IssmDouble>(H_elastic);
|
|---|
| 202 | xDelete<IssmDouble>(H_elastic_local);
|
|---|
| [20029] | 203 | }
|
|---|
| [20036] | 204 |
|
|---|
| [20136] | 205 | /*Transitions: */
|
|---|
| [20690] | 206 | iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
|
|---|
| [20136] | 207 | if(transitions){
|
|---|
| 208 | parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
|
|---|
| 209 |
|
|---|
| 210 | for(int i=0;i<ntransitions;i++){
|
|---|
| 211 | IssmDouble* transition=transitions[i];
|
|---|
| 212 | xDelete<IssmDouble>(transition);
|
|---|
| 213 | }
|
|---|
| 214 | xDelete<IssmDouble*>(transitions);
|
|---|
| [20138] | 215 | xDelete<int>(transitions_M);
|
|---|
| 216 | xDelete<int>(transitions_N);
|
|---|
| [20136] | 217 | }
|
|---|
| 218 |
|
|---|
| [20036] | 219 | /*Requested outputs*/
|
|---|
| [20690] | 220 | iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.slr.requested_outputs");
|
|---|
| [20036] | 221 | if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
|
|---|
| [20690] | 222 | iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
|
|---|
| [20029] | 223 |
|
|---|
| [20036] | 224 |
|
|---|
| [19984] | 225 | }/*}}}*/
|
|---|
| 226 |
|
|---|
| 227 | /*Finite Element Analysis*/
|
|---|
| 228 | void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
|
|---|
| 229 | _error_("not implemented");
|
|---|
| 230 | }/*}}}*/
|
|---|
| 231 | ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
|
|---|
| 232 | /*Default, return NULL*/
|
|---|
| 233 | return NULL;
|
|---|
| 234 | }/*}}}*/
|
|---|
| 235 | ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
|
|---|
| 236 | _error_("Not implemented");
|
|---|
| 237 | }/*}}}*/
|
|---|
| 238 | ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
|
|---|
| 239 | _error_("not implemented yet");
|
|---|
| 240 | }/*}}}*/
|
|---|
| 241 | ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
|
|---|
| 242 | _error_("not implemented yet");
|
|---|
| 243 | }/*}}}*/
|
|---|
| 244 | void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
|
|---|
| 245 | _error_("not implemented yet");
|
|---|
| 246 | }/*}}}*/
|
|---|
| 247 | void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
|
|---|
| 248 | _error_("Not implemented yet");
|
|---|
| 249 | }/*}}}*/
|
|---|
| 250 | void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
|
|---|
| [20153] | 251 |
|
|---|
| 252 | IssmDouble *deltaS = NULL;
|
|---|
| 253 | IssmDouble *S = NULL;
|
|---|
| 254 | int* sidlist = NULL;
|
|---|
| 255 | int numvertices;
|
|---|
| 256 |
|
|---|
| 257 | numvertices= element->GetNumberOfVertices();
|
|---|
| 258 | sidlist=xNew<int>(numvertices);
|
|---|
| 259 |
|
|---|
| 260 | element->GetVerticesSidList(sidlist);
|
|---|
| 261 |
|
|---|
| 262 | deltaS = xNew<IssmDouble>(numvertices);
|
|---|
| 263 | for(int i=0;i<numvertices;i++){
|
|---|
| 264 | deltaS[i]=solution[sidlist[i]];
|
|---|
| 265 | }
|
|---|
| 266 |
|
|---|
| 267 | S = xNew<IssmDouble>(numvertices);
|
|---|
| [20154] | 268 | element->GetInputListOnVertices(S,SealevelEnum,0);
|
|---|
| [20153] | 269 |
|
|---|
| 270 | /*Add deltaS to S:*/
|
|---|
| 271 | for (int i=0;i<numvertices;i++)S[i]+=deltaS[i];
|
|---|
| 272 |
|
|---|
| 273 | /*Add S back into inputs: */
|
|---|
| [20154] | 274 | element->AddInput(SealevelEnum,S,P1Enum);
|
|---|
| [20153] | 275 |
|
|---|
| 276 | /*Free ressources:*/
|
|---|
| 277 | xDelete<int>(sidlist);
|
|---|
| 278 | xDelete<IssmDouble>(deltaS);
|
|---|
| 279 | xDelete<IssmDouble>(S);
|
|---|
| 280 |
|
|---|
| [19984] | 281 | }/*}}}*/
|
|---|
| 282 | void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
|
|---|
| 283 | /*Default, do nothing*/
|
|---|
| 284 | return;
|
|---|
| 285 | }/*}}}*/
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
