![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 25566
- Timestamp:
- 09/11/20 16:37:44 (5 years ago)
- Location:
- issm/branches/trunk-larour-SLPS2020/src/c
- Files:
-
- 3 edited
-
analyses/SealevelriseAnalysis.cpp (modified) (2 diffs)
-
classes/Elements/Tria.cpp (modified) (11 diffs)
-
main/issm_post.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/branches/trunk-larour-SLPS2020/src/c/analyses/SealevelriseAnalysis.cpp
r25535 r25566 175 175 IssmDouble planetradius=0; 176 176 IssmDouble planetarea=0; 177 bool elastic=false; 177 178 178 179 int numoutputs; … … 230 231 } /*}}}*/ 231 232 /*Deal with elasticity {{{*/ 232 /*love numbers: */ 233 iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); 234 iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); 235 iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); 236 iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); 237 iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); 238 iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); 239 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); 240 241 parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); 242 parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); 243 parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); 244 parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); 245 parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); 246 parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); 247 248 /*compute elastic green function for a range of angles*/ 249 iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); 250 M=reCast<int,IssmDouble>(180./degacc+1.); 251 252 // AD performance is sensitive to calls to ensurecontiguous. 253 // // Providing "t" will cause ensurecontiguous to be called. 254 #ifdef _HAVE_AD_ 255 G_rigid=xNew<IssmDouble>(M,"t"); 256 G_elastic=xNew<IssmDouble>(M,"t"); 257 U_elastic=xNew<IssmDouble>(M,"t"); 258 H_elastic=xNew<IssmDouble>(M,"t"); 259 #else 260 G_rigid=xNew<IssmDouble>(M); 261 G_elastic=xNew<IssmDouble>(M); 262 U_elastic=xNew<IssmDouble>(M); 263 H_elastic=xNew<IssmDouble>(M); 264 #endif 265 266 /*compute combined legendre + love number (elastic green function:*/ 267 m=DetermineLocalSize(M,IssmComm::GetComm()); 268 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); 269 // AD performance is sensitive to calls to ensurecontiguous. 270 // // Providing "t" will cause ensurecontiguous to be called. 271 #ifdef _HAVE_AD_ 272 G_elastic_local=xNew<IssmDouble>(m,"t"); 273 G_rigid_local=xNew<IssmDouble>(m,"t"); 274 U_elastic_local=xNew<IssmDouble>(m,"t"); 275 H_elastic_local=xNew<IssmDouble>(m,"t"); 276 #else 277 G_elastic_local=xNew<IssmDouble>(m); 278 G_rigid_local=xNew<IssmDouble>(m); 279 U_elastic_local=xNew<IssmDouble>(m); 280 H_elastic_local=xNew<IssmDouble>(m); 281 #endif 282 283 for(int i=lower_row;i<upper_row;i++){ 284 IssmDouble alpha,x; 285 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0; 286 287 G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); 288 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; 289 U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; 290 H_elastic_local[i-lower_row]= 0; 291 IssmDouble Pn = 0.; 292 IssmDouble Pn1 = 0.; 293 IssmDouble Pn2 = 0.; 294 IssmDouble Pn_p = 0.; 295 IssmDouble Pn_p1 = 0.; 296 IssmDouble Pn_p2 = 0.; 297 298 for (int n=0;n<nl;n++) { 299 IssmDouble deltalove_G; 300 IssmDouble deltalove_U; 301 302 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]); 303 deltalove_U = (love_h[n]-love_h[nl-1]); 304 305 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */ 306 if(n==0){ 307 Pn=1; 308 Pn_p=0; 233 iomodel->FetchData(&elastic,"md.solidearth.settings.elastic"); 234 if(elastic){ 235 /*love numbers: */ 236 iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); 237 iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); 238 iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); 239 iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); 240 iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); 241 iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); 242 parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); 243 244 parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); 245 parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); 246 parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); 247 parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); 248 parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); 249 parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); 250 251 /*compute elastic green function for a range of angles*/ 252 iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); 253 M=reCast<int,IssmDouble>(180./degacc+1.); 254 255 // AD performance is sensitive to calls to ensurecontiguous. 256 // // Providing "t" will cause ensurecontiguous to be called. 257 #ifdef _HAVE_AD_ 258 G_rigid=xNew<IssmDouble>(M,"t"); 259 G_elastic=xNew<IssmDouble>(M,"t"); 260 U_elastic=xNew<IssmDouble>(M,"t"); 261 H_elastic=xNew<IssmDouble>(M,"t"); 262 #else 263 G_rigid=xNew<IssmDouble>(M); 264 G_elastic=xNew<IssmDouble>(M); 265 U_elastic=xNew<IssmDouble>(M); 266 H_elastic=xNew<IssmDouble>(M); 267 #endif 268 269 /*compute combined legendre + love number (elastic green function:*/ 270 m=DetermineLocalSize(M,IssmComm::GetComm()); 271 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); 272 // AD performance is sensitive to calls to ensurecontiguous. 273 // // Providing "t" will cause ensurecontiguous to be called. 274 #ifdef _HAVE_AD_ 275 G_elastic_local=xNew<IssmDouble>(m,"t"); 276 G_rigid_local=xNew<IssmDouble>(m,"t"); 277 U_elastic_local=xNew<IssmDouble>(m,"t"); 278 H_elastic_local=xNew<IssmDouble>(m,"t"); 279 #else 280 G_elastic_local=xNew<IssmDouble>(m); 281 G_rigid_local=xNew<IssmDouble>(m); 282 U_elastic_local=xNew<IssmDouble>(m); 283 H_elastic_local=xNew<IssmDouble>(m); 284 #endif 285 286 for(int i=lower_row;i<upper_row;i++){ 287 IssmDouble alpha,x; 288 alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0; 289 290 G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); 291 G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; 292 U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; 293 H_elastic_local[i-lower_row]= 0; 294 IssmDouble Pn = 0.; 295 IssmDouble Pn1 = 0.; 296 IssmDouble Pn2 = 0.; 297 IssmDouble Pn_p = 0.; 298 IssmDouble Pn_p1 = 0.; 299 IssmDouble Pn_p2 = 0.; 300 301 for (int n=0;n<nl;n++) { 302 IssmDouble deltalove_G; 303 IssmDouble deltalove_U; 304 305 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]); 306 deltalove_U = (love_h[n]-love_h[nl-1]); 307 308 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */ 309 if(n==0){ 310 Pn=1; 311 Pn_p=0; 312 } 313 else if(n==1){ 314 Pn = cos(alpha); 315 Pn_p = 1; 316 } 317 else{ 318 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n; 319 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n; 320 } 321 Pn2=Pn1; Pn1=Pn; 322 Pn_p2=Pn_p1; Pn_p1=Pn_p; 323 324 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential 325 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement 326 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements 309 327 } 310 else if(n==1){ 311 Pn = cos(alpha); 312 Pn_p = 1; 313 } 314 else{ 315 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n; 316 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n; 317 } 318 Pn2=Pn1; Pn1=Pn; 319 Pn_p2=Pn_p1; Pn_p1=Pn_p; 320 321 G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential 322 U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement 323 H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements 324 } 328 } 329 330 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/ 331 int* recvcounts=xNew<int>(IssmComm::GetSize()); 332 int* displs=xNew<int>(IssmComm::GetSize()); 333 334 //recvcounts: 335 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); 336 337 /*displs: */ 338 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); 339 340 /*All gather:*/ 341 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); 342 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); 343 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); 344 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); 345 /*free ressources: */ 346 xDelete<int>(recvcounts); 347 xDelete<int>(displs); 348 349 /*}}}*/ 350 351 /*Avoid singularity at 0: */ 352 G_rigid[0]=G_rigid[1]; 353 parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); 354 G_elastic[0]=G_elastic[1]; 355 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); 356 U_elastic[0]=U_elastic[1]; 357 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); 358 H_elastic[0]=H_elastic[1]; 359 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); 360 361 /*free ressources: */ 362 xDelete<IssmDouble>(love_h); 363 xDelete<IssmDouble>(love_k); 364 xDelete<IssmDouble>(love_l); 365 xDelete<IssmDouble>(love_th); 366 xDelete<IssmDouble>(love_tk); 367 xDelete<IssmDouble>(love_tl); 368 xDelete<IssmDouble>(G_rigid); 369 xDelete<IssmDouble>(G_rigid_local); 370 xDelete<IssmDouble>(G_elastic); 371 xDelete<IssmDouble>(G_elastic_local); 372 xDelete<IssmDouble>(U_elastic); 373 xDelete<IssmDouble>(U_elastic_local); 374 xDelete<IssmDouble>(H_elastic); 375 xDelete<IssmDouble>(H_elastic_local); 325 376 } 326 327 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/328 int* recvcounts=xNew<int>(IssmComm::GetSize());329 int* displs=xNew<int>(IssmComm::GetSize());330 331 //recvcounts:332 ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());333 334 /*displs: */335 ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());336 337 /*All gather:*/338 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());339 ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());340 ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());341 ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());342 /*free ressources: */343 xDelete<int>(recvcounts);344 xDelete<int>(displs);345 346 /*}}}*/347 348 /*Avoid singularity at 0: */349 G_rigid[0]=G_rigid[1];350 parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));351 G_elastic[0]=G_elastic[1];352 parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));353 U_elastic[0]=U_elastic[1];354 parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));355 H_elastic[0]=H_elastic[1];356 parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));357 358 /*free ressources: */359 xDelete<IssmDouble>(love_h);360 xDelete<IssmDouble>(love_k);361 xDelete<IssmDouble>(love_l);362 xDelete<IssmDouble>(love_th);363 xDelete<IssmDouble>(love_tk);364 xDelete<IssmDouble>(love_tl);365 xDelete<IssmDouble>(G_rigid);366 xDelete<IssmDouble>(G_rigid_local);367 xDelete<IssmDouble>(G_elastic);368 xDelete<IssmDouble>(G_elastic_local);369 xDelete<IssmDouble>(U_elastic);370 xDelete<IssmDouble>(U_elastic_local);371 xDelete<IssmDouble>(H_elastic);372 xDelete<IssmDouble>(H_elastic_local);373 377 374 378 /*Indicate we have not yet ran the Geometry Core module: */ -
issm/branches/trunk-larour-SLPS2020/src/c/classes/Elements/Tria.cpp
r25534 r25566 5539 5539 this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); 5540 5540 5541 /*compute area and add to inputs:*/ 5542 area=GetAreaSpherical(); 5543 this->inputs2->SetDoubleInput(AreaEnum,this->lid,area); 5544 this->AddInput2(SealevelAreaEnum,&area,P0Enum); 5545 if(!computerigid & !computeelastic)return; 5546 5541 5547 /*recover precomputed green function kernels:*/ 5542 5548 DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); … … 5561 5567 } 5562 5568 5563 /*compute area:*/ 5564 area=GetAreaSpherical(); 5565 5566 5569 5567 5570 /* Where is the centroid of this element:*/ 5568 5571 … … 5629 5632 this->inputs2->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); 5630 5633 } 5631 this->inputs2->SetDoubleInput(AreaEnum,this->lid,area); 5632 this->AddInput2(SealevelAreaEnum,&area,P0Enum); 5633 5634 5634 5635 /*Free allocations:*/ 5635 5636 #ifdef _HAVE_RESTRICT_ … … 5663 5664 bool notfullygrounded=false; 5664 5665 bool scaleoceanarea= false; 5666 bool computeelastic= true; 5665 5667 5666 5668 /*output: */ … … 5704 5706 #endif 5705 5707 5706 /*recover material parameters:*/5708 /*recover some parameters:*/ 5707 5709 rho_ice=FindParam(MaterialsRhoIceEnum); 5708 5710 rho_water=FindParam(MaterialsRhoSeawaterEnum); 5709 5710 /*recover ocean area scaling: */ 5711 this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); 5711 5712 this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); 5712 5713 5713 5714 /*retrieve precomputed G:*/ 5714 this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);5715 if(computeelastic)this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); 5715 5716 5716 5717 /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ … … 5763 5764 _assert_(!xIsNan<IssmDouble>(bslrice)); 5764 5765 5765 /*convert from m to kg/m^2:*/ 5766 I=I*rho_ice*phi; 5767 5768 /*convolve:*/ 5769 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*I; 5766 if(computeelastic){ 5767 /*convert from m to kg/m^2:*/ 5768 I=I*rho_ice*phi; 5769 5770 /*convolve:*/ 5771 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*I; 5772 } 5770 5773 5771 5774 /*return :*/ … … 5782 5785 bool notfullygrounded=false; 5783 5786 bool scaleoceanarea= false; 5787 bool computeelastic= true; 5784 5788 5785 5789 /*elastic green function:*/ … … 5807 5811 /*If we are here, we are on earth, not on ice: */ 5808 5812 5809 /*recover materialparameters: */5813 /*recover parameters: */ 5810 5814 rho_water=FindParam(MaterialsRhoSeawaterEnum); 5811 5815 rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum); 5812 5813 /*recover ocean area scaling: */ 5816 this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); 5814 5817 this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); 5815 5818 5816 5819 /*retrieve precomputed G:*/ 5817 this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);5820 if(computeelastic)this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); 5818 5821 5819 5822 /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ … … 5831 5834 _assert_(!xIsNan<IssmDouble>(bslrhydro)); 5832 5835 5833 /*convert from m to kg/m^2:*/ 5834 W=W*rho_freshwater*phi; 5835 5836 /*convolve:*/ 5837 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*W; 5836 if(computeelastic){ 5837 /*convert from m to kg/m^2:*/ 5838 W=W*rho_freshwater*phi; 5839 5840 /*convolve:*/ 5841 for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*W; 5842 } 5838 5843 5839 5844 /*output:*/ … … 5848 5853 IssmDouble BP; //change in bottom pressure (Farrel and Clarke, Equ. 4) 5849 5854 IssmDouble constant; 5855 bool computeelastic= true; 5850 5856 5851 5857 /*elastic green function:*/ … … 5854 5860 /*water properties: */ 5855 5861 IssmDouble rho_water; 5862 5863 /*exit now?:*/ 5864 this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); 5865 if(!computeelastic)return; 5856 5866 5857 5867 /*we are here to compute fingerprints originating fromn bottom pressure changes:*/ -
issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp
r25538 r25566 2 2 * \brief: ISSM POST processing main program 3 3 */ 4 5 4 /*includes and prototypes:*/ 6 5 #include "./issm.h" … … 693 692 694 693 /*open file: */ 695 _printf0_(" opening file: " << file < "\n");694 _printf0_(" opening file: " << file << "\n"); 696 695 FILE* fid=fopen(file,"rb"); 697 696 if(fid==NULL)_error_("cound not open file: " << file << "\n");
Note:
See TracChangeset
for help on using the changeset viewer.
