Ignore:
Timestamp:
03/30/21 19:32:13 (4 years ago)
Author:
Eric.Larour
Message:

CHG: hooking the new optimized code and removing the old one. Moving everything
away from FemModel into the sea level core (including Ivins deformation model).

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r26114 r26165  
    47484748/*}}}*/
    47494749#endif
    4750 #ifdef _HAVE_SEALEVELCHANGE_
    4751 void FemModel::SealevelchangeBarystatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition,SealevelMasks* masks) { /*{{{*/
    4752 
    4753         /*serialized vectors:*/
    4754         IssmDouble  bslcice       = 0.;
    4755         IssmDouble  bslcice_cpu   = 0.;
    4756         IssmDouble  bslchydro       = 0.;
    4757         IssmDouble  bslchydro_cpu   = 0.;
    4758         IssmDouble  area      = 0.;
    4759         IssmDouble  oceanarea      = 0.;
    4760         IssmDouble  oceanarea_cpu  = 0.;
    4761         int bp_compute_fingerprints= 0;
    4762         bool isoceantransport=false;
    4763 
    4764         Vector<IssmDouble>* bslcice_partition=NULL;
    4765         IssmDouble* bslcice_partition_serial=NULL;
    4766         IssmDouble* partitionice=NULL;
    4767         int npartice,nel;
    4768 
    4769         Vector<IssmDouble>* bslchydro_partition=NULL;
    4770         IssmDouble* bslchydro_partition_serial=NULL;
    4771         IssmDouble* partitionhydro=NULL;
    4772         bool istws=0;
    4773         int nparthydro;
    4774                
    4775         int npartocean;
    4776         Vector<IssmDouble>* bslcocean_partition=NULL;
    4777         IssmDouble* bslcocean_partition_serial=NULL;
    4778         IssmDouble* partitionocean=NULL;
    4779 
    4780    /*Initialize temporary vector that will be used to sum barystatic components
    4781     * on all local elements, prior to assembly:*/
    4782         int gsize = this->nodes->NumberOfDofs(GsetEnum);
    4783         IssmDouble* RSLgi=xNewZeroInit<IssmDouble>(gsize);
    4784         int* indices=xNew<int>(gsize);
    4785    for(int i=0;i<gsize;i++) indices[i]=i;
    4786 
    4787         /*First, figure out the area of the ocean, which is needed to compute the barystatic component: */
    4788         int i = -1;
    4789         for(Object* & object : this->elements->objects){
    4790                 i +=1;
    4791                 Element* element = xDynamicCast<Element*>(object);
    4792                 element->GetInputValue(&area,AreaEnum);
    4793                 if (masks->isoceanin[i]) oceanarea_cpu += area;
    4794         }
    4795         ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4796         ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4797         _assert_(oceanarea>0.);
    4798 
    4799         /*Initialize partition vectors to retrieve barystatic contributions: */
    4800         this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
    4801         if(npartice){
    4802                 this->parameters->FindParam(&partitionice,&nel,NULL,SolidearthPartitionIceEnum);
    4803                 bslcice_partition= new Vector<IssmDouble>(npartice);
    4804         }
    4805 
    4806         this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
    4807         if(nparthydro){
    4808                 this->parameters->FindParam(&partitionhydro,&nel,NULL,SolidearthPartitionHydroEnum);
    4809                 bslchydro_partition= new Vector<IssmDouble>(nparthydro);
    4810         }
    4811 
    4812         this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum);
    4813         if(npartocean){
    4814                 this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum);
    4815                 bslchydro_partition= new Vector<IssmDouble>(npartocean);
    4816         }
    4817         /*For later:
    4818         npartbarystatic=npartice;
    4819         if(nparthydro>npartbarystatic)npartbarystatic=nparthydro;
    4820         if(npartocean>npartbarystatic)npartbarystatic=npartocean;
    4821         bslc_partition=new Matrix(IssmDouble>(npartbarystatic,3);
    4822 
    4823         bslc_cpu[0]=0; bslc_cpu[1]=0; bslc_cpu[2]=0;
    4824         for(Object* & object : this->elements->objects){
    4825                 Element* element = xDynamicCast<Element*>(object);
    4826                 element->SealevelchangeBarystaticLoads(&bslc_cpu[0], localloads,masks, bslcice_partition,partitionice,oceanarea);
    4827         }
    4828         MPI Bcast localloads -> loads
    4829 
    4830         for(Object* & object : this->elements->objects){
    4831                 Element* element = xDynamicCast<Element*>(object);
    4832                 element->SealevelchangeConvolution(loads);
    4833         }
    4834         */
    4835 
    4836 
    4837 
    4838 
    4839 
    4840         /*Call the barystatic sea level change core for ice : */
    4841         bslcice_cpu=0;
    4842         for(Object* & object : this->elements->objects){
    4843                 Element* element = xDynamicCast<Element*>(object);
    4844                 bslcice_cpu+=element->SealevelchangeBarystaticIce(RSLgi,masks, bslcice_partition,partitionice,oceanarea);
    4845         }
    4846 
    4847         /*Call the barystatic sea level change core for hydro: */
    4848         bslchydro_cpu=0; //make sure to initialize this, so we have a total barystatic contribution computed at 0.
    4849         this->parameters->FindParam(&istws,TransientIshydrologyEnum);
    4850         if(istws){
    4851                 for(int i=0;i<elements->Size();i++){
    4852                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4853                         bslchydro_cpu+=element->SealevelchangeBarystaticHydro(RSLgi,masks, bslchydro_partition,partitionhydro,oceanarea);
    4854                 }
    4855         }
    4856 
    4857         /*Call the barystatic sea level change core for bottom pressures: */
    4858         this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum);
    4859         this->parameters->FindParam(&isoceantransport,TransientIsoceantransportEnum);
    4860         if(bp_compute_fingerprints && isoceantransport){
    4861                 for(int i=0;i<elements->Size();i++){
    4862                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4863                         element->SealevelchangeBarystaticBottomPressure(RSLgi,masks);
    4864                 }
    4865         }
    4866 
    4867         /*Plug values once and assemble: */
    4868         pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL);
    4869         pRSLgi->Assemble();
    4870 
    4871         /*Sum all barystatic components from all cpus:*/
    4872         ISSM_MPI_Reduce (&bslcice_cpu,&bslcice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4873         ISSM_MPI_Bcast(&bslcice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4874         _assert_(!xIsNan<IssmDouble>(bslcice));
    4875 
    4876         ISSM_MPI_Reduce (&bslchydro_cpu,&bslchydro,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4877         ISSM_MPI_Bcast(&bslchydro,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4878         _assert_(!xIsNan<IssmDouble>(bslchydro));
    4879 
    4880         /*Take care of partition vectors:*/
    4881         if(bslcice_partition){
    4882                 bslcice_partition->Assemble();
    4883                 bslcice_partition_serial=bslcice_partition->ToMPISerial();
    4884         }
    4885         if(bslchydro_partition){
    4886                 bslchydro_partition->Assemble();
    4887                 bslchydro_partition_serial=bslchydro_partition->ToMPISerial();
    4888         }
    4889 
    4890 
    4891         /*Free ressources:*/
    4892         xDelete<int>(indices);
    4893         xDelete<IssmDouble>(RSLgi);
    4894         if(bslchydro_partition)delete bslchydro_partition;
    4895         if(bslcice_partition)delete bslcice_partition;
    4896         if(partitionhydro)xDelete<IssmDouble>(partitionhydro);
    4897         if(partitionice)xDelete<IssmDouble>(partitionice);
    4898 
    4899         /*Assign output pointers:*/
    4900         *poceanarea = oceanarea;
    4901         *pbslcice  = bslcice;
    4902         *pbslchydro  = bslchydro;
    4903         *pbslc=bslchydro+bslcice;
    4904         *pbslcice_partition=bslcice_partition_serial;
    4905         *pbslchydro_partition=bslchydro_partition_serial;
    4906 
    4907 }
    4908 /*}}}*/
    4909 void FemModel::SealevelchangeSal(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old,  SealevelMasks* masks, bool verboseconvolution){/*{{{*/
    4910 
    4911         /*serialized vectors:*/
    4912         IssmDouble* RSLg_old=NULL;
    4913 
    4914         IssmDouble* RSLgo  = NULL;
    4915         int* indices = NULL;
    4916         int         gsize;
    4917 
    4918         bool computerigid = true;
    4919         bool computeelastic= true;
    4920 
    4921         /*recover computational flags: */
    4922         this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    4923 
    4924         /*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior
    4925          * to assembly:*/
    4926         gsize = this->nodes->NumberOfDofs(GsetEnum);
    4927         RSLgo=xNewZeroInit<IssmDouble>(gsize);
    4928         indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
    4929 
    4930         /*Serialize vectors from previous iteration:*/
    4931         RSLg_old=pRSLg_old->ToMPISerial();
    4932 
    4933         /*Call the sal sea level change core only if required: */
    4934         if(computerigid){
    4935                 for(Object* & object : this->elements->objects){
    4936                         Element* element = xDynamicCast<Element*>(object);
    4937                         element->SealevelchangeSal(RSLgo,RSLg_old,masks);
    4938                 }
    4939         }
    4940         pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);
    4941         pRSLgo->Assemble();
    4942 
    4943         /*Free ressources:*/
    4944         xDelete<int>(indices);
    4945         xDelete<IssmDouble>(RSLgo);
    4946         xDelete<IssmDouble>(RSLg_old);
    4947 
    4948 }
    4949 /*}}}*/
    4950 void FemModel::SealevelchangeRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz,  SealevelMasks* masks){/*{{{*/
    4951 
    4952         /*serialized vectors:*/
    4953         bool spherical=true;
    4954         IssmDouble* RSLg_old=NULL;
    4955         IssmDouble*     tide_love_h  = NULL;
    4956         IssmDouble*     tide_love_k  = NULL;
    4957         IssmDouble*     load_love_k  = NULL;
    4958         IssmDouble  tide_love_k2secular;
    4959         IssmDouble  moi_e, moi_p, omega, g;
    4960         IssmDouble      m1, m2, m3;
    4961         IssmDouble      lati, longi, radi, value;
    4962         IssmDouble          *latitude    = NULL;
    4963         IssmDouble          *longitude    = NULL;
    4964         IssmDouble          *radius    = NULL;
    4965 
    4966         /*Serialize vectors from previous iteration:*/
    4967         RSLg_old=pRSLg_old->ToMPISerial();
    4968 
    4969         IssmDouble moi_list[3]={0,0,0};
    4970         IssmDouble moi_list_cpu[3]={0,0,0};
    4971         for(Object* & object : this->elements->objects){
    4972                 Element* element = xDynamicCast<Element*>(object);
    4973                 element->SealevelchangeMomentOfInertia(&moi_list[0],RSLg_old,masks );
    4974                 moi_list_cpu[0] += moi_list[0];
    4975                 moi_list_cpu[1] += moi_list[1];
    4976                 moi_list_cpu[2] += moi_list[2];
    4977         }
    4978         ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4979         ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4980         //
    4981         ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4982         ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4983         //
    4984         ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4985         ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4986 
    4987         /*pull out some useful parameters: */
    4988         parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
    4989         parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
    4990         parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
    4991         parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
    4992         parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
    4993         parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
    4994         parameters->FindParam(&omega,RotationalAngularVelocityEnum);
    4995 
    4996         /*compute perturbation terms for angular velocity vector: */
    4997         m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0];
    4998         m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1];
    4999         m3 = -(1+load_love_k[2])/moi_p * moi_list[2];   // term associated with fluid number (3-order-of-magnitude smaller) is negelected
    5000 
    5001         /*recover lat,long and radius vectors from vertices: */
    5002         VertexCoordinatesx(&latitude,&longitude,&radius,this->vertices,spherical);
    5003 
    5004         /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11)
    5005          * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3)
    5006          * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered)
    5007          * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180]
    5008          */
    5009         for(Object* & object : vertices->objects){
    5010                 Vertex* vertex=xDynamicCast<Vertex*>(object);
    5011                 int sid=vertex->Sid();
    5012 
    5013                 lati=latitude[sid]/180*PI;
    5014                 longi=longitude[sid]/180*PI;
    5015                 radi=radius[sid];
    5016 
    5017                 /*only first order terms are considered now: */
    5018                 value=((1.0+tide_love_k[2]-tide_love_h[2])/9.81)*pow(omega*radi,2.0)*
    5019                                                 (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi)));
    5020 
    5021                 pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
    5022         }
    5023 
    5024         /*Assemble mesh velocity*/
    5025         pRSLgo_rot->Assemble();
    5026 
    5027         /*Assign output pointers:*/
    5028         if(pIxz)*pIxz=moi_list[0];
    5029         if(pIyz)*pIyz=moi_list[1];
    5030         if(pIzz)*pIzz=moi_list[2];
    5031         xDelete<IssmDouble>(latitude);
    5032         xDelete<IssmDouble>(longitude);
    5033         xDelete<IssmDouble>(tide_love_h);
    5034         xDelete<IssmDouble>(tide_love_k);
    5035         xDelete<IssmDouble>(load_love_k);
    5036 
    5037         xDelete<IssmDouble>(radius);
    5038 
    5039         /*Free ressources:*/
    5040         xDelete<IssmDouble>(RSLg_old);
    5041 
    5042 }
    5043 /*}}}*/
    5044 void FemModel::SealevelchangeDeformation(Vector<IssmDouble>* pgeoid, Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, SealevelMasks* masks){/*{{{*/
    5045 
    5046         /*serialized vectors:*/
    5047         IssmDouble* RSLg=NULL;
    5048 
    5049         IssmDouble* Up  = NULL;
    5050         IssmDouble* North  = NULL;
    5051         IssmDouble* East  = NULL;
    5052         int* indices = NULL;
    5053         int  gsize;
    5054         int  horiz;
    5055 
    5056         /*retrieve parameters:*/
    5057         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    5058 
    5059         /*Serialize vectors from previous iteration:*/
    5060         RSLg=pRSLg->ToMPISerial();
    5061 
    5062         /*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior
    5063          * to assembly:*/
    5064         gsize = this->nodes->NumberOfDofs(GsetEnum);
    5065         Up=xNewZeroInit<IssmDouble>(gsize);
    5066         if(horiz){
    5067                 North=xNewZeroInit<IssmDouble>(gsize);
    5068                 East=xNewZeroInit<IssmDouble>(gsize);
    5069         }
    5070         indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
    5071 
    5072         /*Call the deformation from loading routines:*/
    5073         for(Object* & object : this->elements->objects){
    5074                 Element* element = xDynamicCast<Element*>(object);
    5075                 element->DeformationFromSurfaceLoads(Up,North,East,RSLg,masks);
    5076         }
    5077 
    5078         pUp->SetValues(gsize,indices,Up,ADD_VAL);
    5079         pUp->Assemble();
    5080 
    5081 
    5082         /*Add RSL to Up to find the geoid: */
    5083         pUp->Copy(pgeoid); pgeoid->AXPY(pRSLg,1);
    5084 
    5085         if (horiz){
    5086                 pNorth->SetValues(gsize,indices,North,ADD_VAL);
    5087                 pNorth->Assemble();
    5088                 pEast->SetValues(gsize,indices,East,ADD_VAL);
    5089                 pEast->Assemble();
    5090         }
    5091 
    5092         /*Free ressources:*/
    5093         xDelete<IssmDouble>(Up);
    5094         if(horiz){
    5095                 xDelete<IssmDouble>(North);
    5096                 xDelete<IssmDouble>(East);
    5097         }
    5098         xDelete<int>(indices);
    5099         xDelete<IssmDouble>(RSLg);
    5100 }
    5101 /*}}}*/
    5102 IssmDouble FemModel::SealevelchangeOceanAverage(Vector<IssmDouble>* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/
    5103 
    5104         IssmDouble* RSLg_serial=NULL;
    5105         IssmDouble  oceanvalue,oceanvalue_cpu;
    5106 
    5107         /*Serialize vectors from previous iteration:*/
    5108         RSLg_serial=RSLg->ToMPISerial();
    5109 
    5110         /*Initialize:*/
    5111         oceanvalue_cpu=0;
    5112 
    5113         /*Go through elements, and add contribution from each element and divide by overall ocean area:*/
    5114         for(Object* & object : this->elements->objects){
    5115                 Element* element = xDynamicCast<Element*>(object);
    5116                 oceanvalue_cpu += element->OceanAverage(RSLg_serial,masks);
    5117         }
    5118 
    5119         ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    5120         ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    5121 
    5122         /*Free ressources:*/
    5123         xDelete<IssmDouble>(RSLg_serial);
    5124 
    5125         return oceanvalue/oceanarea;
    5126 }
    5127 /*}}}*/
    5128 void FemModel::IvinsDeformation(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/
    5129 
    5130         /*Find the litho material to be used by all the elements:*/
    5131         Matlitho* matlitho=NULL;
    5132         for (Object* & object: this->materials->objects){
    5133                 Material* material=xDynamicCast<Material*>(object);
    5134                 if(material->ObjectEnum()==MatlithoEnum){
    5135                         matlitho=xDynamicCast<Matlitho*>(material);
    5136                         break;
    5137                 }
    5138         }
    5139         _assert_(matlitho);
    5140 
    5141 
    5142         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    5143         for(Object* & object : this->elements->objects){
    5144                 Element* element = xDynamicCast<Element*>(object);
    5145                 element->GiaDeflection(wg,dwgdt, matlitho, x,y);
    5146         }
    5147 
    5148         /*Assemble parallel vector:*/
    5149         dwgdt->Assemble();
    5150         wg->Assemble();
    5151 }
    5152 /*}}}*/
    5153 #endif
    51544750void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/
    51554751
Note: See TracChangeset for help on using the changeset viewer.