Index: /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp (revision 27308) @@ -35,5 +35,5 @@ /*Plug inputs into element:*/ - iomodel->FetchDataToInput(inputs,elements,"md.hydrology.spcwatercolum", HydrologyTwsSpcEnum); + iomodel->FetchDataToInput(inputs,elements,"md.hydrology.spcwatercolumn", HydrologyTwsSpcEnum); /*Initialize sea level cumulated sea level loads :*/ @@ -53,5 +53,4 @@ /*Now, do we really want Tws?*/ if(hydrology_model!=HydrologyTwsEnum) return; - parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); Index: /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp (revision 27308) @@ -20,9 +20,25 @@ void LoveAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ IssmDouble* frequencies = NULL; - int nfreq; + IssmDouble* hypergeom_z = NULL; + IssmDouble* hypergeom_table1 = NULL; + IssmDouble* hypergeom_table2 = NULL; + int nfreq,nz, nalpha; iomodel->FetchData(&nfreq,"md.love.nfreq"); + iomodel->FetchData(&nz,"md.love.hypergeom_nz"); + iomodel->FetchData(&nalpha,"md.love.hypergeom_nalpha"); iomodel->FetchData(&frequencies,NULL,NULL,"md.love.frequencies"); + iomodel->FetchData(&hypergeom_z,NULL,NULL,"md.love.hypergeom_z"); + iomodel->FetchData(&hypergeom_table1,NULL,NULL,"md.love.hypergeom_table1"); + iomodel->FetchData(&hypergeom_table2,NULL,NULL,"md.love.hypergeom_table2"); + parameters->AddObject(new DoubleVecParam(LoveFrequenciesEnum,frequencies,nfreq)); + parameters->AddObject(new DoubleVecParam(LoveHypergeomZEnum,hypergeom_z,nz)); + parameters->AddObject(new DoubleMatParam(LoveHypergeomTable1Enum,hypergeom_table1,nz,nalpha)); + parameters->AddObject(new DoubleMatParam(LoveHypergeomTable2Enum,hypergeom_table2,nz,nalpha)); + xDelete(frequencies); + xDelete(hypergeom_z); + xDelete(hypergeom_table1); + xDelete(hypergeom_table2); parameters->AddObject(iomodel->CopyConstantObject("md.love.nfreq",LoveNfreqEnum)); @@ -37,5 +53,7 @@ parameters->AddObject(iomodel->CopyConstantObject("md.love.underflow_tol",LoveUnderflowTolEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.love.pw_threshold",LovePostWidderThresholdEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.love.integration_steps_per_layer",LoveIntStepsPerLayerEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.min_integration_steps",LoveMinIntegrationStepsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.max_integration_dr",LoveMaxIntegrationdrEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.integration_scheme",LoveIntegrationSchemeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.love.istemporal",LoveIsTemporalEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.love.n_temporal_iterations",LoveNTemporalIterationsEnum)); @@ -45,8 +63,12 @@ parameters->AddObject(iomodel->CopyConstantObject("md.love.core_mantle_boundary",LoveCoreMantleBoundaryEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.love.complex_computation",LoveComplexComputationEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.quad_precision",LoveQuadPrecisionEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.debug",LoveDebugEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.hypergeom_nz",LoveHypergeomNZEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.hypergeom_nalpha",LoveHypergeomNAlphaEnum)); }/*}}}*/ void LoveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 27308) @@ -49,4 +49,6 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); + iomodel->FetchDataToInput(inputs,elements,"md.solidearth.transfercount",CouplingTransferCountEnum); + /*external solidearthsolution: solid-Earth model*/ iomodel->FetchData(&isexternal,"md.solidearth.isexternal"); @@ -59,5 +61,5 @@ /*Resolve Mmes using the modelid, if necessary:*/ - if (inputs->GetInputObjectEnum(SolidearthExternalDisplacementEastRateEnum)==DatasetInputEnum){ + if (inputs->GetInputObjectEnum(SolidearthExternalDisplacementUpRateEnum)==DatasetInputEnum){ int modelid; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27307) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27308) @@ -406,5 +406,5 @@ virtual void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0; virtual void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0; - virtual void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids)=0; + virtual void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount)=0; virtual void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; virtual void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 27307) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 27308) @@ -227,5 +227,5 @@ void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");}; + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids,int* vcount){_error_("not implemented yet");}; void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 27307) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 27308) @@ -180,5 +180,5 @@ void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");}; + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids,int* vcount){_error_("not implemented yet");}; void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 27307) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 27308) @@ -187,5 +187,5 @@ void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");}; + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount){_error_("not implemented yet");}; void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 27308) @@ -2198,4 +2198,8 @@ //compute sea level load weights this->GetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset); + + //failsafe for phi so small that GetFractionGeometry returns 0 + if (phi==0) phi=1e-16; + for (int i=0;iGetBarycenterFromLevelset(platbar,plongbar, phi, fraction1, fraction2, late, longe, point1,istrapeze1,planetradius); @@ -6405,5 +6409,5 @@ } /*}}}*/ -void Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){ /*{{{*/ +void Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* n_activevertices){ /*{{{*/ /*Declarations:{{{*/ @@ -6555,11 +6559,15 @@ } - AlphaIndex=xNewZeroInit(3*nel); - if (horiz) AzimuthIndex=xNewZeroInit(3*nel); + AlphaIndex=xNewZeroInit(n_activevertices[this->lid]*nel); + if (horiz) AzimuthIndex=xNewZeroInit(n_activevertices[this->lid]*nel); int intmax=pow(2,16)-1; + int* activevertices=xNew(n_activevertices[this->lid]); + + int av=0; for (int i=0;i<3;i++){ if(lids[this->vertices[i]->lid]==this->lid){ + activevertices[av]=i; for(int e=0;e(intmax*(atan2(dy,dx)/2/M_PI)); + AzimuthIndex[av*nel+e]=reCast(intmax*(atan2(dy,dx)/2/M_PI)); } - AlphaIndex[i*nel+e]=index; - } + AlphaIndex[av*nel+e]=index; + } + av+=1; } //for (int i=0;i<3;i++) } //for(int e=0;einputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*3); - if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*3); + this->inputs->SetIntArrayInput(SealevelchangeConvolutionVerticesEnum,this->lid,activevertices,n_activevertices[this->lid]); + this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*n_activevertices[this->lid]); + if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*n_activevertices[this->lid]); /*}}}*/ @@ -6713,4 +6723,5 @@ /*Free allocations:{{{*/ #ifdef _HAVE_RESTRICT_ + delete activevertices; delete AlphaIndex; if(horiz) AzimuthIndex; @@ -6726,4 +6737,5 @@ #else + xDelete(activevertices); xDelete(AlphaIndex); if(horiz){ @@ -6757,4 +6769,6 @@ IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim; IssmDouble xyz_list[NUMVERTICES][3]; + int* activevertices = NULL; + int n_activevertices, av; #ifdef _HAVE_RESTRICT_ @@ -6831,48 +6845,53 @@ if(horiz) AzimIndex=xNew(SLGEOM_NUMLOADS); + this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); + // 0<=n_activevertices<=3 is the number of vertices this element is in charge of computing fields in during the sea level convolutions + // activevertices contains the vertex indices (1,2 and/or 3) in case debugging is required, they are supposed to appear in the same order as slgeom->lids //Allocate: for(int l=0;lnbar[l]; - AlphaIndex[l]=xNewZeroInit(3*nbar); - if(horiz) AzimIndex[l]=xNewZeroInit(3*nbar); - - - for (int i=0;i<3;i++){ - if(slgeom->lids[this->vertices[i]->lid]==this->lid){ - for(int e=0;elatbarycentre[l][e]; - IssmDouble longe= slgeom->longbarycentre[l][e]; - late=late/180*M_PI; - longe=longe/180*M_PI; - - lati=latitude[i]; - longi=longitude[i]; - + AlphaIndex[l]=xNewZeroInit(n_activevertices*nbar); + if(horiz) AzimIndex[l]=xNewZeroInit(n_activevertices*nbar); + + //av=0; + //for (int i=0;i<3;i++){ + for (int av=0;avlids[this->vertices[i]->lid]==this->lid){ + int i=activevertices[av]; + for(int e=0;elatbarycentre[l][e]; + IssmDouble longe= slgeom->longbarycentre[l][e]; + late=late/180*M_PI; + longe=longe/180*M_PI; + lati=latitude[i]; + longi=longitude[i]; if(horiz){ - /*Compute azimuths*/ + /*Compute azimuths*/ dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi); dy=sin(longe-longi)*cos(late); //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax] - AzimIndex[l][i*nbar+e]=reCast(intmax*(atan2(dy,dx)/2/M_PI)); + AzimIndex[l][av*nbar+e]=reCast(intmax*(atan2(dy,dx)/2/M_PI)); } - /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ - delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda; - alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0))); - doubleindex=alpha/M_PI*reCast(M-1); //maps 0(doubleindex); //truncates doubleindex to integer part - - if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour - if (index==M-1){ //avoids out of bound case - index-=1; - lincoef=1; - } - AlphaIndex[l][i*nbar+e]=index; + /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ + delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda; + alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0))); + doubleindex=alpha/M_PI*reCast(M-1); //maps 0(doubleindex); //truncates doubleindex to integer part + + if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour + if (index==M-1){ //avoids out of bound case + index-=1; + lincoef=1; } - } + AlphaIndex[l][av*nbar+e]=index; + //} + //av+=1; + } + } } @@ -6880,6 +6899,6 @@ /*Save all these arrayins for each element:*/ for (int l=0;linputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*3); - if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*3); + this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*n_activevertices); + if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*n_activevertices); } /*}}}*/ @@ -7227,4 +7246,5 @@ for(int i=0;iLoadWeigths[SLGEOM_OCEAN][i][this->lid]=loadweightsocean[i]; + #ifdef _ISSM_DEBUG_ /*{{{*/ /*Inform mask: */ @@ -7337,16 +7357,19 @@ oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]; } + + oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid]; + oceanaverage*=rho_water*oceanarea; + + /*add ocean average in the global sealevelloads vector:*/ + if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ + int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; + loads->vsubsealevelloads->SetValue(intj,oceanaverage,INS_VAL); + loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL); + } + else loads->vsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL); + #ifdef _ISSM_DEBUG_ this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum); #endif - oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid]; - - /*add ocean average in the global sealevelloads vector:*/ - if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ - int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; - loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water*oceanarea,INS_VAL); - loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL); - } - else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water*oceanarea,INS_VAL); /*add ocean area into a global oceanareas vector:*/ @@ -7367,7 +7390,7 @@ IssmDouble* G=NULL; IssmDouble* Grot=NULL; + IssmDouble* rslfield=NULL; DoubleVecParam* parameter; bool computefuture=false; - int spatial_component=0; bool sal = false; @@ -7388,9 +7411,9 @@ parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL); - this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size); - for (int l=0;linputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size); if (rotation) this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size); - this->SealevelchangeGxL(sealevelpercpu, spatial_component=0,AlphaIndex, AlphaIndexsub, NULL, NULL, G, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false); + rslfield=this->SealevelchangeGxL(G,Grot,loads,polarmotionvector,slgeom,nel,computefuture=false); + this->SealevelchangeCollectGrdfield(sealevelpercpu,rslfield,slgeom,nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false); + } @@ -7406,8 +7429,4 @@ int nel,nbar; bool sal = false; - int* AlphaIndex=NULL; - int* AzimIndex=NULL; - int* AlphaIndexsub[SLGEOM_NUMLOADS]; - int* AzimIndexsub[SLGEOM_NUMLOADS]; int spatial_component=0; IssmDouble* G=NULL; @@ -7418,4 +7437,5 @@ IssmDouble* GNrot=NULL; IssmDouble* GErot=NULL; + IssmDouble* grdfield=NULL; DoubleVecParam* parameter; @@ -7438,7 +7458,4 @@ if(sal){ - this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size); - for (int l=0;linputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size); - parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL); @@ -7449,7 +7466,4 @@ if(horiz){ - this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size); - for (int l=0;linputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size); - parameter = static_cast(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&GH,NULL); @@ -7464,12 +7478,21 @@ } } - - this->SealevelchangeGxL(&RSLGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL,G, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true); + //Relative sea level convolution + grdfield=this->SealevelchangeGxL(G,Grot,loads,polarmotionvector,slgeom,nel,computefuture=true); + this->SealevelchangeCollectGrdfield(&RSLGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true); if(elastic){ - this->SealevelchangeGxL(&UGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL, GU, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true); + //Bedrock Uplift + grdfield=this->SealevelchangeGxL(GU,GUrot,loads,polarmotionvector,slgeom,nel,computefuture=true); + this->SealevelchangeCollectGrdfield(&UGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true); + if(horiz){ - this->SealevelchangeGxL(&NGrd[0],spatial_component=1,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true); - this->SealevelchangeGxL(&EGrd[0],spatial_component=2,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true); + //Bedrock North displacement + grdfield=this->SealevelchangeHorizGxL(spatial_component=1,GH,GNrot,loads,polarmotionvector,slgeom,nel,computefuture=true); + this->SealevelchangeCollectGrdfield(&NGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true); + + //Bedrock East displacement + grdfield=this->SealevelchangeHorizGxL(spatial_component=2,GH,GErot,loads,polarmotionvector,slgeom,nel,computefuture=true); + this->SealevelchangeCollectGrdfield(&EGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true); } } @@ -7496,18 +7519,15 @@ } /*}}}*/ -void Tria::SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/ +IssmDouble* Tria::SealevelchangeGxL(IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture) { /*{{{*/ //This function performs the actual convolution between Green functions and surface Loads for a particular grd field - + int* AlphaIndex=NULL; + int* AlphaIndexsub[SLGEOM_NUMLOADS]; + int* activevertices=NULL; IssmDouble* grdfield=NULL; - int i,e,l,t,a, index, nbar; + int i,e,l,t,a, index, nbar, size, av,ae,b,c; bool rotation=false; - IssmDouble* Centroid_loads=NULL; - IssmDouble* Centroid_loads_copy=NULL; - IssmDouble* Subelement_loads[SLGEOM_NUMLOADS]; - IssmDouble* Subelement_loads_copy[SLGEOM_NUMLOADS]; - IssmDouble* horiz_projection=NULL; - IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS]; int nt=1; //important, ensures there is a defined value if computeviscous is false + int n_activevertices=0; //viscous @@ -7515,4 +7535,240 @@ int viscousindex=0; //important int viscousnumsteps=1; //important + + this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum); + this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); + + //Get green functions indexing & geometry + this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); //the order in which the vertices appear here should be the same as in slgeom->lids + this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size); + for (int l=0;linputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size); + + if(computeviscous){ + this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); + this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum); + if(computefuture) { + nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity + if (nt>viscousnumsteps) nt=viscousnumsteps; + } + else nt=1; + } + //allocate + grdfield=xNewZeroInit(3*nt); + + //early return + if (n_activevertices==0) return grdfield; + + if(rotation){ //add rotational feedback + for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; + b=i*nt; + for (int m=0;m<3;m++){ //polar motion components + for(t=0;tlids[this->vertices[i]->lid]!=this->lid)continue; + b=i*nt; + c=av*nel; + for (int ae=0;aenactiveloads;ae++){ + e=loads->combined_loads_index[ae]; + a=AlphaIndex[c+e]*viscousnumsteps; + for(t=0;tcombined_loads[ae]; + } + } + for(l=0;lnbar[l]; + c=av*nbar; + for (ae=0;aenactivesubloads[l];ae++){ + e=loads->combined_subloads_index[l][ae]; + a=AlphaIndexsub[l][c+e]*viscousnumsteps; + for(t=0;tcombined_subloads[l][ae]; + } + } + } + //av+=1; + } /*}}}*/ + + return grdfield; + +} /*}}}*/ +IssmDouble* Tria::SealevelchangeHorizGxL(int spatial_component, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture) { /*{{{*/ + + //This function performs the actual convolution between Green functions and surface Loads for a particular grd field + int* AlphaIndex=NULL; + int* AzimIndex=NULL; + int* AlphaIndexsub[SLGEOM_NUMLOADS]; + int* AzimIndexsub[SLGEOM_NUMLOADS]; + int* activevertices = NULL; + IssmDouble* grdfield=NULL; + int i,e,l,t,a,b,c, index, nbar, av, ae,n_activevertices, size; + bool rotation=false; + IssmDouble* projected_loads=NULL; + IssmDouble* projected_subloads[SLGEOM_NUMLOADS]; + IssmDouble* horiz_projection=NULL; + IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS]; + int nt=1; //important, ensures there is a defined value if computeviscous is false + + //viscous + bool computeviscous=false; + int viscousindex=0; //important + int viscousnumsteps=1; //important + + //Get green functions indexing & geometry + this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); + this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size); + for (int l=0;linputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size); + this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size); + for (int l=0;linputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size); + + //First, figure out how many time steps to compute grdfield for + this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum); + this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); + if(computeviscous){ + this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); + this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum); + if(computefuture) { + nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity + if (nt>viscousnumsteps) nt=viscousnumsteps; + } + else nt=1; + } + //allocate + grdfield=xNewZeroInit(3*nt); + if (n_activevertices==0) return grdfield; + + if(rotation){ //add rotational feedback + for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; + for (int m=0;m<3;m++){ //polar motion components + for(t=0;t(loads->nactiveloads); + projected_loads=xNewZeroInit(loads->nactiveloads); + for(l=0;lnbar[l]; + projected_subloads[l]=xNewZeroInit(loads->nactivesubloads[l]); + horiz_projectionsub[l]=xNewZeroInit(loads->nactivesubloads[l]); + } + + + //Convolution + //av=0; + for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; + b=i*nt; + + //GxL needs to be projected on the right axis before summation into the grd field + //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency + + //get projection + if (spatial_component==1){ //north + for (int ae=0;aenactiveloads;ae++){ + e=loads->combined_loads_index[ae]; + horiz_projection[ae]=cos(2.0*M_PI*reCast(AzimIndex[av*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int + } + for(l=0;lnbar[l]; + for (ae=0;aenactivesubloads[l];ae++){ + e=loads->combined_subloads_index[l][ae]; + horiz_projectionsub[l][ae]=cos(2.0*M_PI*reCast(AzimIndexsub[l][av*nbar+e])/65535.0); + } + } + } + else if (spatial_component==2){ //east + for (int ae=0;aenactiveloads;ae++){ + e=loads->combined_loads_index[ae]; + horiz_projection[ae]=sin(2.0*M_PI*reCast(AzimIndex[av*nel+e])/65535.0); + } + for(l=0;lnbar[l]; + for (ae=0;aenactivesubloads[l];ae++){ + e=loads->combined_subloads_index[l][ae]; + horiz_projectionsub[l][ae]=sin(2.0*M_PI*reCast(AzimIndexsub[l][av*nbar+e])/65535.0); + } + } + } + + //project load in the right direction + for (int ae=0;aenactiveloads;ae++){ + projected_loads[ae]=loads->combined_loads[ae]*horiz_projection[ae]; + } + for(l=0;lnbar[l]; + for (ae=0;aenactivesubloads[l];ae++){ + projected_subloads[l][ae]=loads->combined_subloads[l][ae]*horiz_projectionsub[l][ae]; + } + } + + //do the convolution + c=av*nel; + for (int ae=0;aenactiveloads;ae++){ + e=loads->combined_loads_index[ae]; + a=AlphaIndex[c+e]*viscousnumsteps; + for(t=0;tnbar[l]; + c=av*nbar; + for (ae=0;aenactivesubloads[l];ae++){ + e=loads->combined_subloads_index[l][ae]; + a=AlphaIndexsub[l][c+e]*viscousnumsteps; + for(t=0;t(horiz_projection); + xDelete(projected_loads); + for(l=0;l(projected_subloads[l]); + xDelete(horiz_projectionsub[l]); + } + return grdfield; + +} /*}}}*/ + + +void Tria::SealevelchangeCollectGrdfield(IssmDouble* grdfieldout, IssmDouble* grdfield, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/ + + //This function aligns grdfield with the requested output format: in a size 3 vector or in a size numberofvertices vector + // if compute viscous is on, we also interpolate the field timewise given the current timestepping as well as collect viscous deformation and update the viscous deformation time series for future time steps + int i,e,l,t,a, index, nbar, av, n_activevertices; + int nt=1; + + + //viscous + bool computeviscous=false; + int viscousindex=0; //important + int viscousnumsteps=1; //important + int* activevertices = NULL; IssmDouble* viscousfield=NULL; IssmDouble* grdfieldinterp=NULL; @@ -7522,6 +7778,8 @@ IssmDouble timeacc; + //parameters & initialization this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum); - this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); + this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); + if(computeviscous){ this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); @@ -7538,200 +7796,88 @@ else nt=1; } - //allocate - grdfield=xNewZeroInit(3*nt); - - if(rotation){ //add rotational feedback - for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ - for (int m=0;m<3;m++){ //polar motion components - for(t=0;t(nel); - for (e=0;eloads[e]; - } - for(l=0;lnbar[l]; - Subelement_loads[l]=xNewZeroInit(nbar); - for (e=0;esubloads[l][e]); - } - } - if(loads->sealevelloads){ - for (e=0;esealevelloads[e]); - } - nbar=slgeom->nbar[SLGEOM_OCEAN]; - for (e=0;esubsealevelloads[e]); - } - } - - //Copy loads if dealing with a horizontal component: the result will need to be projected against the North or East axis for each vertex - if (spatial_component!=0){ - horiz_projection=xNewZeroInit(3*nel); - Centroid_loads_copy=xNewZeroInit(nel); - for (e=0;enbar[l]; - Subelement_loads_copy[l]=xNewZeroInit(nbar); - horiz_projectionsub[l]=xNewZeroInit(3*nbar); - for (e=0;elids[this->vertices[i]->lid]!=this->lid)continue; - - if (spatial_component!=0){ //horizontals /*{{{*/ - //GxL needs to be projected on the right axis before summation into the grd field - //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency - if (spatial_component==1){ //north - for (e=0;e(AzimIndex[i*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int - } - for(l=0;lnbar[l]; - for (e=0;e(AzimIndexsub[l][i*nbar+e])/65535.0);; - } - } - } - else if (spatial_component==2){ //east - for (e=0;e(AzimIndex[i*nel+e])/65535.0); - } - for(l=0;lnbar[l]; - for (e=0;e(AzimIndexsub[l][i*nbar+e])/65535.0);; - } - } - } - for (e=0;enbar[l]; - for (e=0;enbar[l]; - for (e=0;elids[this->vertices[i]->lid]!=this->lid)continue; + grdfieldout[this->vertices[i]->lid]=grdfield[i]; + //} + } + } + else{ + for(i=0;i(grdfield); + return; + } + else { //viscous case // we need to do up to 3 things (* = only if computefuture) - // 1*: add new grdfield contribution to the viscous stack for future time steps - // 2: collect viscous grdfield from past loads due at present-day and add it to grdfield[current_time] + // 1: collect viscous grdfield from past loads due at present-day and add it to grdfield[current_time] + // 2*: add new grdfield contribution to the viscous stack for future time steps // 3*: subtract from viscous stack the grdfield that has already been accounted for so we don't add it again at the next time step + + /*update grdfield at present time using viscous stack at present time: */ + for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; + grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex]; + } + /* Map new grdfield generated by present-day loads onto viscous time vector*/ if(computefuture){ //viscousindex time and first time step of grdfield coincide, so just copy that value - for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; - grdfieldinterp[i*viscousnumsteps+viscousindex]= grdfield[i*nt+0]; + for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; + grdfieldinterp[i*viscousnumsteps+viscousindex]=grdfield[i*nt+0]; } if(viscoustimes[viscousindex]lids[this->vertices[i]->lid]!=this->lid)continue; - grdfieldinterp[i*viscousnumsteps+t] = (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)] - +lincoeff*grdfield[i*nt+(t-viscousindex)]; + for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; + int i_time1= i*nt-viscousindex; + int i_time2= i*viscousnumsteps; + for(int t=viscousindex+1;tlids[this->vertices[i]->lid]!=this->lid)continue; - grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex]; - } - - /*update viscous stack with future deformation from present load: */ - if(computefuture){ - for(int t=viscousnumsteps-1;t>=viscousindex;t--){ //we need to go backwards so as not to zero out viscousfield[i*viscousnumsteps+viscousindex] until the end - for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; - //offset viscousfield to remove all deformation that has already been added - viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t] - -grdfieldinterp[i*viscousnumsteps+viscousindex] - -viscousfield[i*viscousnumsteps+viscousindex]; - } - } /*Save viscous stack now that we updated the values:*/ this->inputs->SetArrayInput(viscousenum,this->lid,viscousfield,3*viscousnumsteps); - - /*Free resources:*/ - xDelete(grdfieldinterp); - } - } - /*}}}*/ - - /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/ - if(percpu){ - for(i=0;ilids[this->vertices[i]->lid]==this->lid){ + } + + /*store values computed on vertices*/ + if(percpu){ + for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; grdfieldout[this->vertices[i]->lid]=grdfield[i*nt+0]; - } - } - } - else{ - for(i=0;i(grdfield); - xDelete(Centroid_loads); - for(l=0;l(Subelement_loads[l]); - if (spatial_component!=0){ - xDelete(horiz_projection); - xDelete(Centroid_loads_copy); - for(l=0;l(Subelement_loads_copy[l]); - xDelete(horiz_projectionsub[l]); - } - } - if (computeviscous){ + //} + } + } + else{ + for(i=0;i(grdfield); xDelete(viscoustimes); if (computefuture){ xDelete(grdfieldinterp); } - } - + /*}}}*/ + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 27307) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 27308) @@ -172,5 +172,5 @@ #ifdef _HAVE_SEALEVELCHANGE_ void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); - void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids); + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount); void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom); void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae); @@ -245,5 +245,7 @@ void UpdateConstraintsExtrudeFromBase(void); void UpdateConstraintsExtrudeFromTop(void); - void SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture); + IssmDouble* SealevelchangeGxL(IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture); + IssmDouble* SealevelchangeHorizGxL(int spatial_component, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture); + void SealevelchangeCollectGrdfield(IssmDouble* grdfieldout, IssmDouble* grdfield, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 27308) @@ -784,4 +784,7 @@ if(hydrology_model==HydrologypismEnum){ analyses_temp[numanalyses++]=HydrologyPismAnalysisEnum; + } + if(hydrology_model==HydrologyTwsEnum){ + analyses_temp[numanalyses++]=HydrologyTwsAnalysisEnum; } if(hydrology_model==HydrologydcEnum){ Index: /issm/trunk-jpl/src/c/classes/GrdLoads.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/GrdLoads.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/classes/GrdLoads.cpp (revision 27308) @@ -20,6 +20,11 @@ /*allocate:*/ + nactiveloads=0; + vloads=new Vector(nel); - for (int i=0;i(slgeom->nbar[i]); + for (int i=0;i(slgeom->nbar[i]); + combined_subloads[i]=xNewZeroInit(slgeom->nbar[i]); + } vsealevelloads=new Vector(nel); @@ -27,4 +32,12 @@ vsubsealevelloads=new Vector(slgeom->nbar[SLGEOM_OCEAN]); + + combined_loads=NULL; + combined_loads_index=NULL; + for (int i=0;i(loads); - for(int i=0;i(subloads[i]); - } delete vsealevelloads; xDelete(sealevelloads); delete vsubsealevelloads; xDelete(subsealevelloads); + if (combined_loads) xDelete(combined_loads); + if (combined_loads_index) xDelete(combined_loads_index); + for(int i=0;i(subloads[i]); + if (combined_subloads[i]) xDelete(combined_subloads[i]); + if (combined_subloads_index[i]) xDelete(combined_subloads_index[i]); + } + + }; /*}}}*/ @@ -139,3 +158,105 @@ ISSM_MPI_Bcast(°2coeff[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); -} /*}}}*/ +}; /*}}}*/ +void GrdLoads::Combineloads(int nel,SealevelGeometry* slgeom){ /*{{{*/ + + int e,l, nbar, ae; + //Determine loads /*{{{*/ + nactiveloads=0; + + if (combined_loads) xDelete(combined_loads); + if (combined_loads_index) xDelete(combined_loads_index); + + //find non zero centroid loads, combine with sealevelloads + if(sealevelloads){ + for (e=0;e(nactiveloads); + combined_loads_index=xNewZeroInit(nactiveloads); + + ae=0; + if(sealevelloads){ + for (e=0;enbar[l]; + if (combined_subloads[l]) xDelete(combined_subloads[l]); + if (combined_subloads_index[l]) xDelete(combined_subloads_index[l]); + + //find non zero subelement loads, combine with subsealevelloads + if(subsealevelloads && l==SLGEOM_OCEAN){ + for (e=0;e(nactivesubloads[l]); + combined_subloads_index[l]=xNewZeroInit(nactivesubloads[l]); + + ae=0; + if(subsealevelloads && l==SLGEOM_OCEAN){ + for (e=0;enbar[l]; + for (e=0;enbar[SLGEOM_OCEAN]; + for (e=0;e* vloads=NULL; IssmDouble* loads=NULL; @@ -22,4 +24,8 @@ Vector* vsubsealevelloads=NULL; IssmDouble* subsealevelloads=NULL; + int* combined_loads_index=NULL; + int* combined_subloads_index[SLGEOM_NUMLOADS]; + IssmDouble* combined_loads=NULL; + IssmDouble* combined_subloads[SLGEOM_NUMLOADS]; GrdLoads(int nel, SealevelGeometry* slgeom); @@ -30,4 +36,5 @@ void BroadcastSealevelLoads(void); void SHDegree2Coefficients(IssmDouble* deg2coeff, FemModel* femmodel, SealevelGeometry* slgeom); + void Combineloads(int nel, SealevelGeometry* slgeom); }; #endif /* _SEALEVELGRDLOADS_H_ */ Index: /issm/trunk-jpl/src/c/classes/Inputs/ArrayInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ArrayInput.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/classes/Inputs/ArrayInput.cpp (revision 27308) @@ -44,8 +44,6 @@ ArrayInput* output = new ArrayInput(this->numberofelements_local); - output->N = xNew(this->numberofelements_local); xMemCpy(output->N,this->N,this->numberofelements_local); - output->values = xNew(this->numberofelements_local); for(int i=0;inumberofelements_local;i++){ if(this->values[i]){ Index: /issm/trunk-jpl/src/c/classes/Inputs/IntArrayInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/IntArrayInput.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/classes/Inputs/IntArrayInput.cpp (revision 27308) @@ -44,8 +44,6 @@ IntArrayInput* output = new IntArrayInput(this->numberofelements_local); - output->N = xNew(this->numberofelements_local); xMemCpy(output->N,this->N,this->numberofelements_local); - output->values = xNew(this->numberofelements_local); for(int i=0;inumberofelements_local;i++){ if(this->values[i]){ Index: /issm/trunk-jpl/src/c/cores/love_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 27308) @@ -12,22 +12,54 @@ #ifdef _HAVE_MPLAPACK_ #include +#include #include "mpblas__Float128.h" #include "mplapack__Float128.h" +#include "mplapack_utils__Float128.h" #endif +template IssmDouble DownCastVarToDouble(doubletype var); // pure declaration + +template <> IssmDouble DownCastVarToDouble(IssmDouble var){ + return var; +} +template <> IssmDouble DownCastVarToDouble(IssmComplex var){ + return std::real(var); +} #ifdef _HAVE_MPLAPACK_ -_Float128 a = 0.2345234534512079875620048770134538q; +template <> IssmDouble DownCastVarToDouble<__float128>(__float128 var){ + return static_cast(var); +} +__float128 pow(__float128 x, int y){ + return powq(x,y); +} +__float128 pow(__float128 x, double y){ + return powq(x,y); +} +__float128 pow(double x, __float128 y){ + return powq(x,y); +} + +ostream& operator<<(ostream& os, __float128 x){ + char buf[128]; + quadmath_snprintf(buf, sizeof(buf), "%.34Qf", x); + + os << buf; + return os; +} #endif /*local definitions:*/ -class LoveVariables{ /*{{{*/ +template class LoveVariables{ /*{{{*/ public: - IssmDouble g0; - IssmDouble r0; - IssmDouble* EarthMass; - int nyi; + doubletype g0; + doubletype r0; + doubletype* EarthMass; + int nyi, ifreq, nfreq; int starting_layer; int* deg_layer_delete; + int* nstep; + doubletype* mu; + doubletype* la; LoveVariables(){ /*{{{*/ @@ -35,9 +67,14 @@ r0=0; EarthMass=NULL; + mu=NULL; + la=NULL; nyi=0; + nfreq=0; + ifreq=0; starting_layer=0; + nstep=NULL; deg_layer_delete=NULL; } /*}}}*/ - LoveVariables(IssmDouble* EarthMassin,IssmDouble g0in,IssmDouble r0in,int nyiin,int starting_layerin, int* deg_layer_deletein){ /*{{{*/ + LoveVariables(doubletype* EarthMassin,doubletype g0in,doubletype r0in,int nyiin,int starting_layerin, int* deg_layer_deletein, int* nstepin){ /*{{{*/ EarthMass=EarthMassin; g0=g0in; @@ -46,6 +83,16 @@ starting_layer=starting_layerin; deg_layer_delete=deg_layer_deletein; + nstep=nstepin; + mu=NULL; + la=NULL; + nfreq=0; + ifreq=0; } /*}}}*/ - ~LoveVariables(){}; + ~LoveVariables(){ + xDelete(deg_layer_delete); + xDelete(nstep); + if(mu) xDelete(mu); + if(la) xDelete(la); + }; }; /*}}}*/ @@ -57,5 +104,5 @@ doubletype* L; doubletype* Kernels; - int sh_nmin, sh_nmax, nfreq, nkernels; + int sh_nmin, sh_nmax, nfreq, nkernels, lower_row, nfreqtotal; LoveNumbers(){ /*{{{*/ @@ -69,8 +116,10 @@ nkernels=0; } /*}}}*/ - LoveNumbers(int sh_nminin, int sh_nmaxin, int nfreqin, Matlitho* matlitho){ /*{{{*/ + LoveNumbers(int sh_nminin, int sh_nmaxin, int nfreqin, int lower_rowin,int nfreqtotalin,Matlitho* matlitho){ /*{{{*/ sh_nmax=sh_nmaxin; sh_nmin=sh_nminin; nfreq=nfreqin; + lower_row=lower_rowin; + nfreqtotal=nfreqtotalin; nkernels=(sh_nmax+1)*(matlitho->numlayers+1)*6; H=xNewZeroInit(nfreq*(sh_nmax+1)); @@ -81,8 +130,12 @@ void DownCastToDouble(LoveNumbers* LoveDouble){ for(int i=0;i<(sh_nmax+1)*nfreq;i++){ - LoveDouble->H[i]=std::real(H[i]); - LoveDouble->K[i]=std::real(K[i]); - LoveDouble->L[i]=std::real(L[i]); - } + LoveDouble->H[i]=DownCastVarToDouble(H[i]); + LoveDouble->K[i]=DownCastVarToDouble(K[i]); + LoveDouble->L[i]=DownCastVarToDouble(L[i]); + } + for(int i=0;iKernels[i]=DownCastVarToDouble(Kernels[i]); + } + } void LoveMPI_Gather(LoveNumbers* Love_local, int lower_row){ @@ -112,4 +165,42 @@ xDelete(displs); } + void Broadcast(void){ + //Intended only for IssmDouble type + Vector* vH; + Vector* vK; + Vector* vL; + + vH= new Vector((sh_nmax+1)*nfreqtotal); + for(int i=0;i<(sh_nmax+1)*nfreq;i++) vH->SetValue(lower_row*(sh_nmax+1)+i,DownCastVarToDouble(H[i]),INS_VAL); + xDelete(H); + vH->Assemble(); + H=vH->ToMPISerial(); + delete vH; + + vK= new Vector((sh_nmax+1)*nfreqtotal); + for(int i=0;i<(sh_nmax+1)*nfreq;i++) vK->SetValue(lower_row*(sh_nmax+1)+i,DownCastVarToDouble(K[i]),INS_VAL); + xDelete(K); + vK->Assemble(); + K=vK->ToMPISerial(); + delete vK; + + vL= new Vector((sh_nmax+1)*nfreqtotal); + for(int i=0;i<(sh_nmax+1)*nfreq;i++) vL->SetValue(lower_row*(sh_nmax+1)+i,DownCastVarToDouble(L[i]),INS_VAL); + xDelete(L); + vL->Assemble(); + L=vL->ToMPISerial(); + delete vL; + } + void KernelBroadcast(void){ + Vector* vKernels; + vKernels= new Vector(nkernels*nfreqtotal); + for(int i=0;iSetValue(lower_row*nkernels+i,DownCastVarToDouble(Kernels[i]),INS_VAL); + } + xDelete(Kernels); + vKernels->Assemble(); + Kernels=vKernels->ToMPISerial(); + delete vKernels; + } void Copy(LoveNumbers* LoveDup){ for(int i=0;i<(sh_nmax+1)*nfreq;i++){ @@ -128,5 +219,7 @@ xDelete(Kernels); }; -}; /*}}}*/ +}; + + /*}}}*/ /*self contained support routines used by cores below:*/ @@ -139,4 +232,9 @@ return value; } /*}}}*/ +#ifdef _HAVE_MPLAPACK_ +template <> __float128 angular_frequency<__float128>(IssmDouble frequency){ /*{{{*/ + return 2.0*PI*frequency; +} /*}}}*/ +#endif template void allgesv(int* pnyi, int* pnrhs, doubletype* yilocal, int* plda, int* ipiv, doubletype* rhslocal, int* pldb, int* pinfo); template <> void allgesv(int* pnyi, int* pnrhs, IssmDouble* yilocal, int* plda, int* ipiv, IssmDouble* rhslocal, int* pldb, int* pinfo){ /*{{{*/ @@ -147,4 +245,22 @@ //zgesv_(pnyi, pnrhs, yilocal, plda, ipiv, rhslocal, pldb, pinfo); } /*}}}*/ +#ifdef _HAVE_MPLAPACK_ +template <> void allgesv<__float128>(int* pnyi, int* pnrhs, __float128* yilocal, int* plda, int* ipiv, __float128* rhslocal, int* pldb, int* pinfo){ /*{{{*/ + mplapackint nyi=*pnyi; + mplapackint nrhs=*pnrhs; + mplapackint lda=*plda; + mplapackint* qipiv=NULL; + mplapackint ldb=*pldb; + mplapackint info=0; + qipiv=xNewZeroInit(*pnyi); + + Rgesv(nyi, nrhs, yilocal, lda, qipiv, rhslocal, ldb, info); + + for (int i;i=0;i<*pnyi) ipiv[i]=qipiv[i]; + *pinfo=info; + xDelete(qipiv); +} /*}}}*/ +#endif + template doubletype factorial(int n){ /*{{{*/ doubletype prod=1; @@ -219,5 +335,5 @@ int indxi, indf; - IssmDouble PW_test; + doubletype PW_test; IssmDouble PW_threshold; femmodel->parameters->FindParam(&PW_threshold,LovePostWidderThresholdEnum); @@ -261,51 +377,225 @@ Lovet[t*(sh_nmax+1)+d]=LoveM[NTit-1]; }/*}}}*/ -template void GetEarthRheology(doubletype* pla, doubletype* pmu, int layer_index, doubletype omega, Matlitho* matlitho){ /*{{{*/ + +template doubletype HypergeomTableLookup(doubletype z1, doubletype alpha, IssmDouble* h1, IssmDouble* z, int nz, int nalpha){/*{{{*/ + int iz1, iz2, ialpha; + doubletype lincoef; + doubletype hf,h00,h10, h01, h11, za, zd, ha, hb,hc,hd, m0,m1,t; + doubletype dalpha=1.0/(nalpha-1); // alpha table resolution given 0 <= alpha <= 1 + ialpha= static_cast(DownCastVarToDouble(alpha/dalpha)); + lincoef=alpha/dalpha-ialpha;//linear fraction in [0;1] for alpha interpolation + iz1=nz; + for (int i=0;iabs(z1)) { + iz1=i-1; + break; + } + } + + if (iz1<0){ + //1-hf for very small abs(z) tends to 0, and is very log-linear with respect to log(z), so we can simply extrapolate the value of hf via the loglog slope + hf=(1.0-lincoef)*h1[ialpha*nz+0]+lincoef*h1[(ialpha+1)*nz+0]; + hf=1.0- (1.0-hf)*pow(10.0,(log10(abs(z1))-log10(abs(z[0])))); + //hf[0]=1.0; + } + else if (iz1==nz){ + //hf for very large abs(z) tends to 0, and is very log-linear with respect to log(z), so we can simply extrapolate the value of hf via the loglog slope + hf=(1.0-lincoef)*h1[ialpha*nz+nz-1]+lincoef*h1[(ialpha+1)*nz+nz-1]; + hf=hf *pow(10.0,-(log10(abs(z1))-log10(abs(z[nz-1])))); + //hf[0]=0; + } + else{ //cubic spline interpolation + //edge cases: extrapolate 1 point + if (iz1==0){ + za=2.0*z[0]-z[1]; + ha=(1.0-lincoef)*h1[ialpha*nz+0]+lincoef*h1[(ialpha+1)*nz+0]; + ha=1.0- (1.0-ha) *pow(10.0,log10(abs(za))-log10(abs(z[0]))); + } + else { + za=z[iz1-1]; + ha=(1.0-lincoef)*h1[ialpha*nz+iz1-1] + lincoef*h1[(ialpha+1)*nz+iz1-1]; + } + + if (iz1==nz-2){ + zd=2.0*z[nz-1]-z[nz-2]; + hd=(1.0-lincoef)*h1[ialpha*nz+nz-1]+lincoef*h1[(ialpha+1)*nz+nz-1]; + hd=hd *pow(10.0,-(log10(abs(zd))-log10(abs(z[nz-1])))); + } + else { + zd=z[iz1+2]; + hd=(1.0-lincoef)*h1[ialpha*nz+iz1+2]+lincoef*h1[(ialpha+1)*nz+iz1+2]; + } + + hb=(1.0-lincoef)*h1[ialpha*nz+iz1] +lincoef*h1[(ialpha+1)*nz+iz1]; + hc=(1.0-lincoef)*h1[ialpha*nz+iz1+1] +lincoef*h1[(ialpha+1)*nz+iz1+1]; + + //left derivative + m0= 0.5*(z[iz1+1]-z[iz1])*((hc-hb)/(z[iz1+1]-z[iz1]) + (hb-ha)/(z[iz1]-za)); + //right derivative + m1= 0.5*(z[iz1+1]-z[iz1])*((hd-hc)/(zd-z[iz1+1]) + (hc-hb)/(z[iz1+1]-z[iz1])); + + //interpolation abscissa + t=(z1-z[iz1])/(z[iz1+1]-z[iz1]); + + //cubic spline functions + h00=2*pow(t,3)-3*pow(t,2)+1; + h10=pow(t,3)-2*pow(t,2)+t; + h01=-2*pow(t,3)+3*pow(t,2); + h11=pow(t,3)-pow(t,2); + + hf=h00*hb + h10*m0 + h01*hc + h11*m1; + } + return hf; + +}/*}}}*/ + +template doubletype muEBM(int layer_index, doubletype omega, Matlitho* matlitho, FemModel* femmodel); //pure declaration +template <> IssmComplex muEBM(int layer_index, IssmComplex omega, Matlitho* matlitho, FemModel* femmodel){/*{{{*/ + // Initialization + int nz, nalpha, dummy1, dummy2; + IssmComplex mu; + IssmDouble* z=NULL; + IssmDouble* h1=NULL; + IssmDouble* h2=NULL; + IssmComplex hf11, hf12, hf21, hf22; + IssmDouble factor=0; + IssmComplex z1, z2; + IssmComplex U1, U2; + IssmComplex j=1i; + //Matlitho parameters + IssmDouble alpha=matlitho->ebm_alpha[layer_index]; + IssmDouble delta=matlitho->ebm_delta[layer_index]; + IssmDouble taul=matlitho->ebm_taul[layer_index]; + IssmDouble tauh=matlitho->ebm_tauh[layer_index]; + IssmDouble vi=matlitho->viscosity[layer_index]; + IssmDouble mu0=matlitho->lame_mu[layer_index]; + //fetch hypergeometric function tables and parameters + femmodel->parameters->FindParam(&nz,LoveHypergeomNZEnum); + femmodel->parameters->FindParam(&nalpha,LoveHypergeomNAlphaEnum); + femmodel->parameters->FindParam(&z,&dummy1,LoveHypergeomZEnum); + femmodel->parameters->FindParam(&h1,&dummy1,&dummy2,LoveHypergeomTable1Enum); + femmodel->parameters->FindParam(&h2,&dummy1,&dummy2,LoveHypergeomTable2Enum); + omega=omega/j; + + z1= -pow(omega*tauh,2.0); + z2= -pow(omega*taul,2.0); + //Table1 h1 should be 2F1([1 1+alpha], [2+alpha/2], z) + //Table2 h2 should be 2F1([1 0.5+alpha], [1.5+alpha/2], z) + hf11=HypergeomTableLookup(z1, alpha, h1, z, nz, nalpha); + hf21=HypergeomTableLookup(z1, alpha, h2, z, nz, nalpha); + hf12=HypergeomTableLookup(z2, alpha, h1, z, nz, nalpha); + hf22=HypergeomTableLookup(z2, alpha, h2, z, nz, nalpha); + + + //Ivins et al (2020) p11 + U1=(pow(tauh,alpha)-pow(taul,alpha))/alpha-pow(omega,2.0)/(2.0+alpha)*(pow(tauh,2.0+alpha)*hf11-pow(taul,2.0+alpha)*hf12); + U2=(pow(tauh,1.0+alpha)*hf21-pow(taul,1.0+alpha)*hf22)/(1.0+alpha); + + factor= alpha*delta/(pow(tauh,alpha)-pow(taul,alpha)); + U1=(1.0+factor) *U1; + U2=factor*omega*U2 +mu0/vi/omega; + mu=mu0*(U1+j*U2)/(pow(U1,2.0)+pow(U2,2.0)); + omega=omega*j; + + xDelete(z); + xDelete(h1); + xDelete(h2); + return mu; +}/*}}}*/ + +template <> IssmDouble muEBM(int layer_index, IssmDouble omega, Matlitho* matlitho, FemModel* femmodel){/*{{{*/ + // Initialization + int nz, nalpha, dummy1, dummy2; + IssmDouble mu; + IssmDouble* z=NULL; + IssmDouble* h1=NULL; + IssmDouble hf11, hf12; + IssmDouble factor, B, D, z1, z2; + //Matlitho parameters + IssmDouble alpha=matlitho->ebm_alpha[layer_index]; + IssmDouble delta=matlitho->ebm_delta[layer_index]; + IssmDouble taul=matlitho->ebm_taul[layer_index]; + IssmDouble tauh=matlitho->ebm_tauh[layer_index]; + IssmDouble vi=matlitho->viscosity[layer_index]; + IssmDouble mu0=matlitho->lame_mu[layer_index]; + //fetch hypergeometric function tables and parameters + femmodel->parameters->FindParam(&nz,LoveHypergeomNZEnum); + femmodel->parameters->FindParam(&nalpha,LoveHypergeomNAlphaEnum); + femmodel->parameters->FindParam(&z,&dummy1,LoveHypergeomZEnum); + femmodel->parameters->FindParam(&h1,&dummy1,&dummy2,LoveHypergeomTable1Enum); + + z1=-omega*tauh; + z2=-omega*taul; + //Table1 h1 should be 2F1([1 1+alpha], [2+alpha], z) + hf11=HypergeomTableLookup(z1, alpha, h1, z, nz, nalpha); + hf12=HypergeomTableLookup(z2, alpha, h1, z, nz, nalpha); + + //Ivins et al. (2022) p1979 + factor= alpha*delta/(pow(tauh,alpha)-pow(taul,alpha)); + B= factor/(1.0+alpha) *mu0/vi * (pow(tauh,1.0+alpha)*hf11 - pow(taul,1.0+alpha)*hf12); + D= omega*vi/mu0* 1.0/(1.0+omega*vi/mu0*(1.0+delta) -pow(omega*vi/mu0,2.0)*B); + + xDelete(z); + xDelete(h1); + return mu=mu0*D; +}/*}}}*/ +#ifdef _HAVE_MPLAPACK_ +template <> __float128 muEBM<__float128>(int layer_index, __float128 omega, Matlitho* matlitho, FemModel* femmodel){/*{{{*/ + // Initialization + int nz, nalpha, dummy1, dummy2; + IssmDouble* z=NULL; + IssmDouble* h1=NULL; + __float128 mu; + __float128 hf11, hf12; + __float128 factor, B, D, z1, z2; + //Matlitho parameters + __float128 alpha=matlitho->ebm_alpha[layer_index]; + __float128 delta=matlitho->ebm_delta[layer_index]; + __float128 taul=matlitho->ebm_taul[layer_index]; + __float128 tauh=matlitho->ebm_tauh[layer_index]; + __float128 vi=matlitho->viscosity[layer_index]; + __float128 mu0=matlitho->lame_mu[layer_index]; + //fetch hypergeometric function tables and parameters + femmodel->parameters->FindParam(&nz,LoveHypergeomNZEnum); + femmodel->parameters->FindParam(&nalpha,LoveHypergeomNAlphaEnum); + femmodel->parameters->FindParam(&z,&dummy1,LoveHypergeomZEnum); + femmodel->parameters->FindParam(&h1,&dummy1,&dummy2,LoveHypergeomTable1Enum); + + z1=-(omega*tauh); + z2=-(omega*taul); + //Table1 h1 should be 2F1([1 1+alpha], [2+alpha], z) + hf11=HypergeomTableLookup<__float128>(z1, alpha, h1, z, nz, nalpha); + hf12=HypergeomTableLookup<__float128>(z2, alpha, h1, z, nz, nalpha); + + //Ivins et al. (2022) p1979 + //Note: therein, mu(s') = s'*mu~(s'); s'=omega*tauM=omega*vi/mu0 + factor= alpha*delta/(pow(tauh,alpha)-pow(taul,alpha)); + B= factor/(1.0q+alpha) *mu0/vi * (pow(tauh,1.0q+alpha)*hf11 - pow(taul,1.0q+alpha)*hf12); + D= omega*vi/mu0* 1.0q/(1.0q+omega*vi/mu0*(1.0q+delta) -pow(omega*vi/mu0,2.0q)*B); + + xDelete(z); + xDelete(h1); + return mu=mu0*D; +}/*}}}*/ +#endif +template void GetEarthRheology(doubletype* pla, doubletype* pmu, int layer_index, doubletype omega, Matlitho* matlitho, FemModel* femmodel){ /*{{{*/ //returns lame parameters (material rigity) lambda and mu for the right frequency and layer doubletype mu,la; - IssmDouble vi=matlitho->viscosity[layer_index]; - IssmDouble mu0=matlitho->lame_mu[layer_index]; - IssmDouble la0=matlitho->lame_lambda[layer_index]; + doubletype vi=matlitho->viscosity[layer_index]; + doubletype mu0=matlitho->lame_mu[layer_index]; + doubletype la0=matlitho->lame_lambda[layer_index]; int rheo=matlitho->rheologymodel[layer_index]; if (vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid - IssmDouble ka=la0 + 2.0/3.0*mu0; //Bulk modulus + doubletype ka=la0 + 2.0/3.0*mu0; //Bulk modulus if (rheo==2){//EBM - IssmDouble alpha=matlitho->ebm_alpha[layer_index]; - IssmDouble delta=matlitho->ebm_delta[layer_index]; - IssmDouble taul=matlitho->ebm_taul[layer_index]; - IssmDouble tauh=matlitho->ebm_tauh[layer_index]; - IssmDouble hf1,hf2; - doubletype U1,U2; - IssmDouble* a=xNew(2); - IssmDouble* b=xNew(1); - - a[0]=1.0;a[1]=1.0+alpha/2.0; - b[0]=2.0+alpha/2.0; - //hf1=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*tauh,2.0), 0, 0, 15); - //hf2=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*taul,2.0), 0, 0, 15); - //hf1=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*tauh,2.0)); - //hf2=hypergeometric_pFq_({a[0], a[1]}, {b[0]},-pow(omega*taul,2.0)); - U1=(pow(tauh,alpha)-pow(taul,alpha))/alpha-pow(omega,2.0)/(2.0+alpha)*(pow(tauh,2.0+alpha)*hf1-pow(taul,2.0+alpha)*hf2); - - a[0]=1.0;a[1]=.5+alpha/2.0; - b[0]=1.5+alpha/2.0; - //hf1=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*tauh,2.0), 0, 0, 15); - //hf2=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*taul,2.0), 0, 0, 15); - //hf1=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*tauh,2.0)); - //hf2=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*taul,2.0)); - U2=(pow(tauh,1.0+alpha)*hf1-pow(taul,1.0+alpha)*hf2)/(1.0+alpha); - - mu=mu/(1.0+ alpha*delta/(pow(tauh,alpha)-pow(taul,alpha))*(U1 + omega*U2) -mu/vi/omega); + mu=muEBM(layer_index, omega, matlitho, femmodel); la=ka-2.0/3.0*mu; - - xDelete(a); - xDelete(b); } else if (rheo==1){//Burgers - IssmDouble vi2=matlitho->burgers_viscosity[layer_index]; - IssmDouble mu2=matlitho->burgers_mu[layer_index]; + doubletype vi2=matlitho->burgers_viscosity[layer_index]; + doubletype mu2=matlitho->burgers_mu[layer_index]; mu=mu0*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu0/vi)+mu0/vi2*omega); @@ -321,18 +611,40 @@ mu=mu0; } + *pla=la; *pmu=mu; } /*}}}*/ -IssmDouble GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/ + +template void EarthRheology(LoveVariables* vars, IssmDouble* frequencies, int nfreq, Matlitho* matlitho, FemModel* femmodel){/*{{{*/ + doubletype omega; + //reset pointers to NULL if this function was previously called + if(vars->mu) xDelete(vars->mu); + if(vars->la) xDelete(vars->la); + //precompute rheology at the requested frequencies + vars->mu=xNewZeroInit(nfreq*matlitho->numlayers); + vars->la=xNewZeroInit(nfreq*matlitho->numlayers); + vars->nfreq=nfreq; + for (int i=0;inumlayers;i++){ + for (int fr=0;fr(frequencies[fr]); + GetEarthRheology(&vars->la[i*nfreq+fr], &vars->mu[i*nfreq+fr], i,omega,matlitho, femmodel); + //cout << i << " " << fr << " " << vars->mu[i*nfreq+fr] << "\n"; + } + } +}/*}}}*/ + +template doubletype GetGravity(doubletype r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/ //computes gravity at radius r2 - IssmDouble* EarthMass; - IssmDouble g, GG; + doubletype* EarthMass; + doubletype g, GG; + IssmDouble GGp; EarthMass=vars->EarthMass; - femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); - IssmDouble ro=matlitho->density[layer_index]; - IssmDouble M=0; - IssmDouble r1=0; + femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum); + GG=GGp; + doubletype ro=matlitho->density[layer_index]; + doubletype M=0; + doubletype r1=0; if (layer_index==0){ M=4.0/3.0*PI*ro*pow(r2,3.0); @@ -344,21 +656,23 @@ return g= GG*M/pow(r2,2.0); }/*}}}*/ -template void fill_yi_prefactor(doubletype* yi_prefactor, int* pdeg, doubletype* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ +template void fill_yi_prefactor(doubletype* yi_prefactor, int* pdeg, doubletype* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ //precalculates partial derivative factors for function yi_derivatives - IssmDouble ra=matlitho->radius[matlitho->numlayers]; - IssmDouble g0,r0,mu0, GG; - int nstep,nindex, starting_layer; - - femmodel->parameters->FindParam(&mu0,LoveMu0Enum); + doubletype ra=matlitho->radius[matlitho->numlayers]; + doubletype g0,r0,mu0; + IssmDouble mu0p, GG; + int nstep,nsteps,nindex, starting_layer; + + femmodel->parameters->FindParam(&mu0p,LoveMu0Enum); femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); g0=vars->g0; r0=vars->r0; + mu0=mu0p; starting_layer=vars->starting_layer; doubletype frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm; - IssmDouble xmin,xmax,x,dr; - IssmDouble g,ro, issolid; + doubletype xmin,xmax,x,dr; + doubletype g,ro, issolid; if (pomega) { //frequency and degree dependent terms /*{{{*/ @@ -370,9 +684,14 @@ for (int layer_index=starting_layer;layer_indexnumlayers;layer_index++){ + nstep=vars->nstep[layer_index]; + nsteps=0; + for (int i=0;instep[i]; ro=matlitho->density[layer_index]; issolid=matlitho->issolid[layer_index]; if(issolid==1){ - GetEarthRheology(&la, &mu,layer_index,omega,matlitho); + //GetEarthRheology(&la, &mu,layer_index,omega,matlitho); + mu=vars->mu[layer_index*vars->nfreq+vars->ifreq]; + la=vars->la[layer_index*vars->nfreq+vars->ifreq]; /*_______Expressions*/ @@ -386,5 +705,5 @@ f[2]=(la*fn/flm); f[3]=rm0*rlm; - f[4]=ro*pow(omega,2.0)*ra/mu0; + f[4]=-ro*pow(omega,2.0)*ra*ra/mu0; f[5]=(-4.0*mu/flm); f[6]=fn*frh; @@ -399,8 +718,12 @@ dr = (xmax -xmin)/nstep; x=xmin; + + //fixme + g=GetGravity((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars); + for (int n=0;n(x*ra,layer_index,femmodel,matlitho,vars); + nindex=nsteps*36+n*36; yi_prefactor[nindex+ 0*6+0]= f[0]/x; // in dy[0*6+0] yi_prefactor[nindex+ 0*6+1]= f[1]; // in dy[0*6+1] @@ -423,4 +746,8 @@ for (int layer_index=starting_layer;layer_indexnumlayers;layer_index++){ + nstep=vars->nstep[layer_index]; + nsteps=0; + for (int i=0;instep[i]; + ro=matlitho->density[layer_index]; issolid=matlitho->issolid[layer_index]; @@ -433,7 +760,12 @@ dr = (xmax -xmin)/nstep; x=xmin; + + + //fixme + g=GetGravity((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars); + for (int n=0;n(x*ra,layer_index,femmodel,matlitho,vars); if(issolid==1){ @@ -450,4 +782,8 @@ } else { // static terms /*{{{*/ for (int layer_index=starting_layer;layer_indexnumlayers;layer_index++){ + nstep=vars->nstep[layer_index]; + nsteps=0; + for (int i=0;instep[i]; + ro=matlitho->density[layer_index]; issolid=matlitho->issolid[layer_index]; @@ -461,7 +797,9 @@ dr = (xmax -xmin)/nstep; x=xmin; + //fixme + g=GetGravity((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars); for (int n=0;n(x*ra,layer_index,femmodel,matlitho,vars); + nindex=nsteps*36+n*36; if(issolid==1){ yi_prefactor[nindex+ 1*6+5]= -frhg0; // in dy[1*6+5] @@ -484,10 +822,14 @@ } }/*}}}*/ -template void yi_derivatives(doubletype* dydx, doubletype* y, int layer_index, int n, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ +template void yi_derivatives(doubletype* dydx, doubletype* y, int layer_index, int n, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ //computes yi derivatives at r=radius[layer_index]+ n/nstep*(radius[layer_index+1]-radius[layer_index]) - IssmDouble issolid=matlitho->issolid[layer_index]; - int iy,id,ny, nindex, nstep; - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + int issolid=matlitho->issolid[layer_index]; + int iy,id,ny, nindex, nstep, nsteps; + //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + + nstep=vars->nstep[layer_index]; + nsteps=0; + for (int i=0;instep[i]; /*{{{*/ /* For reference: @@ -552,6 +894,8 @@ dy[1*6+0]= (-4.0*(fgr/rg0)+fn/x)/x; dy[1*6+1]= -2.0/x-fgr/rg0; + } */ /*}}}*/ + nindex=nsteps*36+n*36; if(issolid==1){ @@ -564,5 +908,4 @@ dydx[id]=0.0; for (iy=0;iy void propagate_yi_euler(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ - //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer? +template void propagate_yi_euler(doubletype* y, doubletype xmin, doubletype xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ + //computes this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer? //euler method int nstep; - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + nstep=vars->nstep[layer_index]; + doubletype* dydx=xNewZeroInit(6); - IssmDouble dr = (xmax -xmin)/nstep; - IssmDouble x=xmin; + doubletype dr = (xmax -xmin)/nstep; + doubletype x=xmin; for(int i = 0;i(dydx,y,layer_index, i,yi_prefactor,femmodel,matlitho); + yi_derivatives(dydx,y,layer_index, i,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++){ y[j]+=dydx[j]*dr; @@ -588,9 +933,10 @@ xDelete(dydx); }/*}}}*/ -template void propagate_yi_RK2(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ - //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer? +template void propagate_yi_RK2(doubletype* y, doubletype xmin, doubletype xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ + //computes this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer? //Implements Runge-Kutta 2nd order (midpoint method) int nstep; - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + nstep=vars->nstep[layer_index]; doubletype k1[6]={0}; @@ -602,11 +948,11 @@ doubletype y3[6]={0}; - IssmDouble dr = (xmax -xmin)/nstep; - IssmDouble x=xmin; + doubletype dr = (xmax -xmin)/nstep; + doubletype x=xmin; for(int i = 0;i(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho); + yi_derivatives(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;} - yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++){ @@ -616,28 +962,29 @@ } }/*}}}*/ -template void propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ - //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer? + template void propagate_yi_RK4(doubletype* y, doubletype xmin, doubletype xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/ + //computes this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer? //Implements Runge-Kutta 4th order int nstep; - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); - - IssmDouble* k1=xNewZeroInit(6); - IssmDouble* k2=xNewZeroInit(6); - IssmDouble* k3=xNewZeroInit(6); - IssmDouble* k4=xNewZeroInit(6); - IssmDouble* y1=xNewZeroInit(6); - IssmDouble* y2=xNewZeroInit(6); - IssmDouble* y3=xNewZeroInit(6); - - IssmDouble dr = (xmax -xmin)/nstep; - IssmDouble x=xmin; + //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + nstep=vars->nstep[layer_index]; + + doubletype k1[6]={0}; + doubletype k2[6]={0}; + doubletype k3[6]={0}; + doubletype k4[6]={0}; + doubletype y1[6]={0}; + doubletype y2[6]={0}; + doubletype y3[6]={0}; + + doubletype dr = (xmax -xmin)/nstep; + doubletype x=xmin; for(int i = 0;i(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho); + yi_derivatives(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;} - yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;} - yi_derivatives(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + yi_derivatives(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;} - yi_derivatives(k4,y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho); + yi_derivatives(k4,y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++){ @@ -649,11 +996,11 @@ //Last step: we don't know the derivative at xmax, so we will assume the values at xmax-dr int i=nstep/2; - yi_derivatives(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho); + yi_derivatives(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;} - yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;} - yi_derivatives(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + yi_derivatives(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;} - yi_derivatives(k4,y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + yi_derivatives(k4,y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars); for (int j=0;j<6;j++){ @@ -663,41 +1010,49 @@ x = x + 2.0*dr; - xDelete(k1); - xDelete(k2); - xDelete(k3); - xDelete(k4); - xDelete(y1); - xDelete(y2); - xDelete(y3); }/*}}}*/ -template void Innersphere_boundaryconditions(doubletype* yi, int layer_index, int deg, doubletype omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ +template void Innersphere_boundaryconditions(doubletype* yi, int layer_index, int deg, doubletype omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ //fills the boundary conditions at the bottom of layer[layer_index] in yi[0:2][0:5] - IssmDouble r = matlitho->radius[layer_index]; - IssmDouble ra=matlitho->radius[matlitho->numlayers]; - IssmDouble g0,r0,mu0, GG; int nyi; - - femmodel->parameters->FindParam(&mu0,LoveMu0Enum); - femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); + doubletype r = matlitho->radius[layer_index]; + doubletype ra=matlitho->radius[matlitho->numlayers]; + doubletype g0,r0,mu0, GG; + IssmDouble mu0p, GGp; + + + femmodel->parameters->FindParam(&mu0p,LoveMu0Enum); + femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum); g0=vars->g0; r0=vars->r0; + mu0=mu0p; + GG=GGp; nyi=vars->nyi; - IssmDouble ro=matlitho->density[layer_index]; - IssmDouble issolid=matlitho->issolid[layer_index]; - IssmDouble g=GetGravity(r,layer_index,femmodel,matlitho,vars); - doubletype la,mu; - - if (layer_index==0){ - // radius[0] cannot be 0 for numerical reasons, but below our first interface at radius[0] would in reality be the same material as in the first layer - GetEarthRheology(&la, &mu,layer_index,omega,matlitho); - } else { - GetEarthRheology(&la, &mu,layer_index-1,omega,matlitho); - } - - IssmDouble cst = 4.0*PI*GG*ro; - IssmDouble r2=pow(r,2.0); + + doubletype g=GetGravity(r,layer_index,femmodel,matlitho,vars); + doubletype la,mu,ro; + + int i=layer_index-1; + if (layer_index==0) i=layer_index; + + ro=matlitho->density[i]; + + //elastic values + la=matlitho->lame_lambda[i]; + mu=matlitho->lame_mu[i]; + doubletype Kappa=(la+2.0/3.0*mu); + + //update to viscoelastic values + mu=vars->mu[i*vars->nfreq+vars->ifreq]; + la = Kappa-2.0/3.0*mu; + + doubletype cst = 4.0*PI*GG*ro; + doubletype r2=pow(r,2.0); + + //Greff-Lefftz and Legros 1997, p701, analytical solution for incompressible elastic layer for y3, y4, y5 ensuring they =0 at r=0 + //These equations are then divided by r^n for numerical stability at higher degrees + + //all terms in y1 y2 y6 are 0 in that layer because they are of the type r^l with l<0 and would diverge at the origin, that's why we only need 3 equations yi[0+nyi*0]=1.0*r/ra; @@ -725,9 +1080,103 @@ yi[5+nyi*2]=deg/(r*g0); + + /*doubletype vp2 = (la + 2.0*mu)/ro; + yi[0+nyi*0]=1.0*r/ra; + yi[0+nyi*1]=1.0/(r*ra); + yi[0+nyi*2]=0.0; + + yi[1+nyi*0]=(2.0*mu*(deg-1.0-3.0/deg) + cst/3.0*ro*r2)/mu0; + yi[1+nyi*1]=(2.0*mu*(deg-1.0)/r2 + cst/3.0*ro)/mu0; + yi[1+nyi*2]=-ro/mu0; + + yi[2+nyi*0]=(deg+3.0)/(deg*(deg+1.0))*r/ra; + yi[2+nyi*1]=1.0/(deg*r*ra); + yi[2+nyi*2]=0.0; + + yi[3+nyi*0]=2.0*mu*(deg+2.0)/((deg+1.0)*mu0); + yi[3+nyi*1]=2.0*mu*(deg-1.0)/(deg*r2*mu0); + yi[3+nyi*2]=0.0; + + yi[4+nyi*0]=0.0; + yi[4+nyi*1]=0.0; + yi[4+nyi*2]=1.0/(g0*ra); + + yi[5+nyi*0]=-cst*r/g0; + yi[5+nyi*1]=-cst/(r*g0); + yi[5+nyi*2]=deg/(r*g0);*/ + + + }/*}}}*/ -template void build_yi_system(doubletype* yi, int deg, doubletype omega, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars) { /*{{{*/ - - IssmDouble g0,r0,mu0,x,ro1, GG; - int nyi,starting_layer, nstep; +template void Coremantle_boundaryconditions(doubletype* yi, int layer_index, int deg, doubletype omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ + //fills the boundary conditions at the bottom of layer[layer_index] in yi[0:2][0:5] + + int nyi; + doubletype r = matlitho->radius[layer_index]; + doubletype ra=matlitho->radius[matlitho->numlayers]; + doubletype g0,r0,mu0, GG; + IssmDouble mu0p, GGp; + + + femmodel->parameters->FindParam(&mu0p,LoveMu0Enum); + femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum); + + g0=vars->g0; + r0=vars->r0; + mu0=mu0p; + GG=GGp; + nyi=vars->nyi; + + doubletype ro=matlitho->density[layer_index-1]; + + + if (!matlitho->issolid[layer_index]) _error_("Love core error: CMB conditions requested but layer " << layer_index << " (mantle side) is not solid"); + if (matlitho->issolid[layer_index-1]) _error_("Love core error: CMB conditions requested but layer " << layer_index-1 << " (outer core) is solid"); + + doubletype cst = 4.0/3.0*PI*GG*ro; + doubletype rl1=pow(r,deg-1); + + yi[0+nyi*0]=-rl1/cst/ra; + yi[0+nyi*1]=1.0/ra; + yi[0+nyi*2]=0.0; + + yi[1+nyi*0]=0.0; + yi[1+nyi*1]=ro*cst*r/mu0; + yi[1+nyi*2]=0.0; + + yi[2+nyi*0]=0.0; + yi[2+nyi*1]=0.0; + yi[2+nyi*2]=1.0/ra; + + yi[3+nyi*0]=0.0; + yi[3+nyi*1]=0.0; + yi[3+nyi*2]=0.0; + + yi[4+nyi*0]=r*rl1/(g0*ra); + yi[4+nyi*1]=0.0; + yi[4+nyi*2]=0.0; + + yi[5+nyi*0]=2.0*(deg-1)*rl1/g0; + yi[5+nyi*1]=3.0*cst/g0; + yi[5+nyi*2]=0.0; + +}/*}}}*/ +template void build_yi_system(doubletype* yi, int deg, doubletype omega, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars) { /*{{{*/ + + doubletype g0,r0,mu0,x,ro1, GG; + int nyi,starting_layer, nstep; + doubletype xmin,xmax,one,ro,g, ra; + IssmDouble mu0p, GGp; + bool debug; + int ny,is,ii,jj; + int scheme; + int ici = 0; // Index of current interface + int cmb=0; + + femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum); + femmodel->parameters->FindParam(&mu0p,LoveMu0Enum); + femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum); + femmodel->parameters->FindParam(&debug,LoveDebugEnum); + femmodel->parameters->FindParam(&scheme,LoveIntegrationSchemeEnum); g0=vars->g0; @@ -735,11 +1184,7 @@ nyi=vars->nyi; starting_layer=vars->starting_layer; - - femmodel->parameters->FindParam(&mu0,LoveMu0Enum); - femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); - - IssmDouble xmin,xmax,one,ro,g; - IssmDouble ra=matlitho->radius[matlitho->numlayers]; + mu0=mu0p; + GG=GGp; + ra=matlitho->radius[matlitho->numlayers]; for (int i=0;i<6*(matlitho->numlayers+1);i++){ @@ -749,11 +1194,11 @@ } - int ny,is,ii,jj; doubletype ystart[6]; for (int k=0;k<6;k++) ystart[k]=0.0; - int ici = 0; // Index of current interface + + for (int i = starting_layer; inumlayers;i++){ - + ici=i-starting_layer; xmin=matlitho->radius[i]/ra; xmax=(matlitho->radius[i+1])/ra; @@ -774,7 +1219,11 @@ // Numerical Integration - //propagate_yi_euler(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho); - propagate_yi_RK2(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho); - //propagate_yi_RK4(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho); + if (debug) propagate_yi_euler(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho, vars); + else { + if (scheme==0) propagate_yi_euler(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho, vars); + else if (scheme==1) propagate_yi_RK2(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho, vars); + else if (scheme==2) propagate_yi_RK4(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho, vars); + else _error_("Love core error: integration scheme not found"); + } // Boundary Condition matrix - propagation part ii = 6*(ici+1)+is; @@ -794,22 +1243,24 @@ } else { // Boundary Condition matrix - liquid regions ro1=matlitho->density[i]; - g=GetGravity(matlitho->radius[i], i, femmodel,matlitho,vars); + g=GetGravity(matlitho->radius[i], i, femmodel,matlitho,vars); ii = 6*ici; - yi[ii+nyi*(ii+3)] = -1.0; - yi[ii+nyi*(ii+4+3)] = -g0/g; - yi[(ii+1)+nyi*(ii+3)]=-ro1*g*ra/mu0; - yi[(ii+2)+nyi*(ii+1+3)]=-1.0; - yi[(ii+5)+nyi*(ii+3)]= 4.0*PI*GG*ro1*ra/g0; - yi[(ii+4)+nyi*(ii+4+3)]=-1.0; - yi[(ii+5)+nyi*(ii+5+3)]=-1.0; - g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho,vars); + jj = 6*ici+3; + yi[ii+nyi*(jj)] = -1.0; + yi[ii+nyi*(jj+4)] = -g0/g; + yi[(ii+1)+nyi*(jj)]=-ro1*g*ra/mu0; + yi[(ii+2)+nyi*(jj+1)]=-1.0; + yi[(ii+5)+nyi*(jj)]= 4.0*PI*GG*ro1*ra/g0; + yi[(ii+4)+nyi*(jj+4)]=-1.0; + yi[(ii+5)+nyi*(jj+5)]=-1.0; + g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho,vars); ii = 6*(ici+1); - yi[ii+nyi*(ii-1)]=-1.0; - yi[ii+nyi*(ii+1)]=yi[(ii+4)+nyi*(ii+1)]*g0/g; // yi(17,14) solution integration 1 of z5 CMB - yi[ii+nyi*(ii+2)]=yi[(ii+4)+nyi*(ii+2)]*g0/g; // yi(17,15) solution integration 2 of z5 CMB + + yi[ii+nyi*(jj+2)]=-1.0; + yi[ii+nyi*(jj+4)]=yi[(ii+4)+nyi*(jj+4)]*g0/g; // yi(17,14) solution integration 1 of z5 CMB + yi[ii+nyi*(jj+5)]=yi[(ii+4)+nyi*(jj+5)]*g0/g; // yi(17,15) solution integration 2 of z5 CMB // yi(13,..) y1 CMB - yi[(ii+1)+nyi*(ii-1)]=-ro1*g*ra/mu0; - yi[(ii+2)+nyi*(ii)]=-1.0; - yi[(ii+5)+nyi*(ii-1)]= 4.0*PI*GG*ro1*ra/g0; + yi[(ii+1)+nyi*(jj+2)]=-ro1*g*ra/mu0; + yi[(ii+2)+nyi*(jj+3)]=-1.0; + yi[(ii+5)+nyi*(jj+2)]= 4.0*PI*GG*ro1*ra/g0; } ici = ici+1; @@ -817,4 +1268,7 @@ //-- Internal sphere: integration starts here rather than r=0 for numerical reasons + /*if (starting_layer==cmb) Coremantle_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho,vars); + else Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho,vars);*/ + Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho,vars); @@ -834,9 +1288,10 @@ }/*}}}*/ -template void yi_boundary_conditions(doubletype* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type){ /*{{{*/ - - IssmDouble g0,r0,mu0,ra,rb,rc; +template void yi_boundary_conditions(doubletype* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type){ /*{{{*/ + + doubletype g0,r0,mu0,ra,rb,rc; int nyi,icb,cmb,starting_layer; - IssmDouble* EarthMass; + doubletype* EarthMass; + IssmDouble mu0p; g0=vars->g0; @@ -846,8 +1301,9 @@ EarthMass=vars->EarthMass; - femmodel->parameters->FindParam(&mu0,LoveMu0Enum); + femmodel->parameters->FindParam(&mu0p,LoveMu0Enum); femmodel->parameters->FindParam(&icb,LoveInnerCoreBoundaryEnum); femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum); + mu0=mu0p; // In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces ra=matlitho->radius[matlitho->numlayers]; @@ -861,5 +1317,5 @@ } - IssmDouble ro_mean=EarthMass[matlitho->numlayers-1]/(4.0/3.0*PI*pow(ra,3.0)); + doubletype ro_mean=EarthMass[matlitho->numlayers-1]/(4.0/3.0*PI*pow(ra,3.0)); for (int i=0;i<(matlitho->numlayers+1)*6;i++) yi_righthandside[i]=0.0; @@ -917,11 +1373,12 @@ } }/*}}}*/ -template void solve_yi_system(doubletype* loveh, doubletype* lovel, doubletype* lovek, int deg, doubletype omega, IssmDouble* frequencies, doubletype* yi, doubletype* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu){ /*{{{*/ - - IssmDouble g0,r0,mu0,loveratio,underflow_tol; +template void solve_yi_system(doubletype* loveh, doubletype* lovel, doubletype* lovek, int deg, doubletype omega, IssmDouble* frequencies, doubletype* yi, doubletype* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu){ /*{{{*/ + + doubletype g0,r0,mu0; //IssmDouble* frequencies; - int nyi,starting_layer, dummy; - bool allow_layer_deletion; - IssmDouble* EarthMass=NULL; + int nyi,starting_layer, dummy,cmb; + bool allow_layer_deletion, debug; + doubletype* EarthMass=NULL; + IssmDouble mu0p,loveratio,underflow_tol; g0=vars->g0; @@ -931,9 +1388,12 @@ EarthMass=vars->EarthMass; - femmodel->parameters->FindParam(&mu0,LoveMu0Enum); + femmodel->parameters->FindParam(&mu0p,LoveMu0Enum); femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum); femmodel->parameters->FindParam(&underflow_tol,LoveUnderflowTolEnum); + femmodel->parameters->FindParam(&debug,LoveDebugEnum); + femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum); //femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); - IssmDouble ra=matlitho->radius[matlitho->numlayers]; + mu0=mu0p; + doubletype ra=matlitho->radius[matlitho->numlayers]; bool exit=false; int lda,ldb; @@ -953,4 +1413,20 @@ } } + + if (debug){ + IssmDouble* yidebug=xNew(nyi*nyi); + IssmDouble* rhsdebug=xNew(nyi); + for (int i=0;iresults->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveYiEnum,yidebug,nyi,nyi,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveRhsEnum,rhsdebug,nyi,1,0,0)); + xDelete(yidebug); + xDelete(rhsdebug); + } + //-- Resolution int* ipiv=xNewZeroInit(nyi); //pivot index vector @@ -974,6 +1450,6 @@ _error_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");*/ - if(VerboseSolution() && info!=0){ - _printf_("i j yi[i+nyi*j] rhs[i]"); + if(VerboseSolution() && verbosecpu && info!=0){ + _printf_("i j yi[i+nyi*j] rhs[i]\n"); for (int i=0;i(frequencies[0]) || deg==0){ @@ -1022,5 +1499,5 @@ if (deg==vars->deg_layer_delete[starting_layer]){ // if we are at the degree where we should delete the current layer, proceed to delete the bottom layer - if (omega!=0 && VerboseSolution() && verbosecpu) _printf_(", deleting layer " << starting_layer << "\n"); + //if (omega!=0 && VerboseSolution() && verbosecpu) _printf_(", deleting layer " << starting_layer << "\n"); nyi-=6; starting_layer+=1; @@ -1035,5 +1512,8 @@ } - Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho,vars); //we move the first interface to the new starting layer. yi[0:2,0:5] will be different + /*if (starting_layer==cmb) Coremantle_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho,vars); + else Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho,vars); //we move the first interface to the new starting layer. yi[0:2,0:5] will be different + */ + Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho,vars); } else { //we are ready to save the outputs and break the main loop @@ -1134,15 +1614,15 @@ postwidder_transform(pmtf_orthot,pmtf_orthof,2,t,2,NTit,xi,femmodel); } + xDelete(pmtf_colinearf); + xDelete(pmtf_orthof); if(VerboseSolution() && verbosecpu) _printf_("done!\n"); } xDelete(xi); - xDelete(pmtf_colinearf); - xDelete(pmtf_orthof); }/*}}}*/ -template void compute_love_numbers(LoveNumbers* Lovef, LoveNumbers* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu){ - - int nstep, kernel_index,kernel_indexe,deleted_layer_offset, deg, sh_nmin, sh_nmax, nfreq; +template void compute_love_numbers(LoveNumbers* Lovef, LoveNumbers* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu){ + + int nsteps, kernel_index,kernel_indexe,deleted_layer_offset, deg, sh_nmin, sh_nmax, nfreq; doubletype lovek, loveh, lovel, loveratio; doubletype omega; @@ -1150,9 +1630,13 @@ doubletype* yi_righthandside=NULL; doubletype* yi=NULL; - IssmDouble underflow_tol; + doubletype underflow_tol; + IssmDouble dr; bool freq_skip, istemporal; - - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + int cmb=0; + int nyi_init=0; + + //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum); + femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum); nfreq=Lovef->nfreq; @@ -1162,17 +1646,21 @@ // reset deleted layers in case we have called this function before; - vars->starting_layer=0; - vars->nyi=6*(matlitho->numlayers+1); - - yi_prefactor=xNewZeroInit(6*6*nstep*matlitho->numlayers); - yi_righthandside=xNewZeroInit(vars->nyi); - yi=xNewZeroInit(vars->nyi*vars->nyi); - - //precalculate yi coefficients that do not depend on degree or frequency - fill_yi_prefactor(yi_prefactor, NULL, NULL,femmodel, matlitho,vars); + vars->starting_layer=0; + vars->nyi=6*(matlitho->numlayers-vars->starting_layer+1); + nyi_init=6*(matlitho->numlayers+1); + nsteps=0; + for (int i=0;inumlayers;i++) nsteps+=vars->nstep[i]; + + //yi_prefactor=xNewZeroInit(6*6*nstep*matlitho->numlayers); + yi_prefactor=xNewZeroInit(6*6*nsteps); + yi_righthandside=xNewZeroInit(nyi_init); + yi=xNewZeroInit(nyi_init*nyi_init); + + //precompute yi coefficients that do not depend on degree or frequency + fill_yi_prefactor(yi_prefactor, NULL, NULL,femmodel, matlitho,vars); if (VerboseSolution() && Elastic && verbosecpu) _printf_("\n"); - for(int deg=0;deg<2;deg++){ // calculation is in the center of mass reference frame, neutralize degree 0 and 1 mass changes + for(int deg=0;deg<2;deg++){ // calculation is in the center of mass reference frame, neutralize degree 0 and 1 mass changes, i.e 1+k=0 for (int fr=0;frK[fr*(sh_nmax+1)+deg]=-1.0; @@ -1185,17 +1673,19 @@ } - //precalculate yi coefficients that depend on degree but not frequency + //precompute yi coefficients that depend on degree but not frequency fill_yi_prefactor(yi_prefactor, °, NULL,femmodel, matlitho,vars); for (int fr=0;fr(frequencies[fr]); + vars->ifreq=fr; - //precalculate yi coefficients that depend on degree and frequency + //precompute yi coefficients that depend on degree and frequency fill_yi_prefactor(yi_prefactor, °,&omega,femmodel, matlitho,vars); //solve the system - yi_boundary_conditions(yi_righthandside,deg,femmodel,matlitho,vars,forcing_type); + yi_boundary_conditions(yi_righthandside,deg,femmodel,matlitho,vars,forcing_type); build_yi_system(yi,deg,omega,yi_prefactor,femmodel,matlitho,vars); - solve_yi_system(&loveh,&lovel,&lovek, deg, omega, frequencies, yi, yi_righthandside,femmodel, matlitho,vars,verbosecpu); + solve_yi_system(&loveh,&lovel,&lovek, deg, omega, frequencies, yi, yi_righthandside,femmodel, matlitho,vars,verbosecpu && !Elastic); + Lovef->H[fr*(sh_nmax+1)+deg]=loveh; Lovef->K[fr*(sh_nmax+1)+deg]=lovek-1.0; @@ -1209,5 +1699,5 @@ } - if (Elastic) { // if elastic values were provided, we copy elastic love numbers above the cutoff degree instead of calculating them + if (Elastic) { // if elastic values were provided, we copy elastic love numbers above the cutoff degree instead of computing them for(int deg=sh_cutoff+1;deg(yi); @@ -1236,5 +1728,5 @@ /*templated cores:*/ -LoveVariables* love_init(FemModel* femmodel, Matlitho* matlitho){/*{{{*/ +template LoveVariables* love_init(FemModel* femmodel, Matlitho* matlitho, bool verbosecpu){/*{{{*/ /*initialize Planet_Mass(r) for efficient computation of gravity, value of surface gravity and inital size of the yi equation system*/ @@ -1242,16 +1734,36 @@ bool verbosemod = (int)VerboseModule(); int numlayers = matlitho->numlayers; - IssmDouble* r = matlitho->radius; - IssmDouble r1,r2,ro, GG; + int minsteps; + doubletype* r=NULL; + doubletype r1,r2,ro, GG; + IssmDouble GGp; + IssmDouble dr; /*outputs:*/ - IssmDouble* EarthMass=NULL; - IssmDouble g0,r0; - int nyi,starting_layer; + doubletype* EarthMass=NULL; + doubletype g0,r0; + int nyi,starting_layer,cmb; int* deg_layer_delete; - - femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); - EarthMass=xNewZeroInit(numlayers+1); + int* nstep; + + + femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum); + femmodel->parameters->FindParam(&minsteps, LoveMinIntegrationStepsEnum); + femmodel->parameters->FindParam(&dr, LoveMaxIntegrationdrEnum); + femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum); + GG=GGp; + EarthMass=xNewZeroInit(numlayers+1); deg_layer_delete=xNewZeroInit(numlayers); + + r=xNewZeroInit(numlayers+1); + nstep=xNewZeroInit(numlayers); + for (int i=0;iradius[i]; + if (iradius[i+1]-matlitho->radius[i])/dr/2)*2; + if (nstep[i](r); + return new LoveVariables(EarthMass,g0,r0,nyi,starting_layer,deg_layer_delete,nstep); } /*}}}*/ @@ -1290,4 +1803,5 @@ bool save_results; bool complex_computation; + bool quad_precision; bool verbosecpu=false; @@ -1295,7 +1809,4 @@ doubletype lovek, loveh, lovel, loveratio; IssmDouble pw_threshold, pw_test_h, pw_test_l,pw_test_k; - - LoveNumbers* Lovef=NULL; - LoveNumbers* Tidalf=NULL; /* parallel computing */ @@ -1309,5 +1820,5 @@ IssmDouble* frequencies_fluid=NULL; - LoveVariables* vars=NULL; + LoveVariables* vars=NULL; /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/ @@ -1332,16 +1843,10 @@ femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum); femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum); + femmodel->parameters->FindParam(&quad_precision,LoveQuadPrecisionEnum); femmodel->parameters->FindParam(&pw_threshold,LovePostWidderThresholdEnum); if (istemporal) femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum); - /*Initialize three love matrices: geoid, vertical and horizontal displacement*/ - Lovef= new LoveNumbers(sh_nmin,sh_nmax,nfreq, matlitho); - Tidalf= new LoveNumbers(2,2,nfreq, matlitho); - Elastic= new LoveNumbers(sh_nmin,sh_nmax,1,matlitho); - Fluid= new LoveNumbers(sh_nmin,sh_nmax,1,matlitho); - - /*Initialize love kernels (real and imaginary parts): */ - vars=love_init(femmodel,matlitho); - + Elastic= new LoveNumbers(sh_nmin,sh_nmax,1,1,1,matlitho); + Fluid= new LoveNumbers(sh_nmin,sh_nmax,1,1,1,matlitho); //distribute frequencies for parallel computation /*{{{*/ int nt_local, nf_local, lower_row, upper_row; @@ -1364,10 +1869,11 @@ if (lower_row==0) verbosecpu=true; //let only cpu1 be verbose if(VerboseSolution() && verbosecpu) _printf0_(" computing LOVE numbers\n"); + vars=love_init(femmodel,matlitho,verbosecpu); frequencies_local=xNewZeroInit(nf_local); for (int fr=0;fr(sh_nmin,sh_nmax,nf_local, matlitho); - Tidalf_local= new LoveNumbers(2,2,nf_local, matlitho); + Lovef_local= new LoveNumbers(sh_nmin,sh_nmax,nf_local,lower_row,nfreq, matlitho); + Tidalf_local= new LoveNumbers(2,2,nf_local,lower_row,nfreq, matlitho); /*}}}*/ @@ -1381,7 +1887,9 @@ // run elastic and fluid love numbers if(VerboseSolution() && verbosecpu) _printf_(" elastic\n"); + EarthRheology(vars,frequencies_elastic,1,matlitho,femmodel); compute_love_numbers(Elastic, NULL, forcing_type, sh_nmax,frequencies_elastic, femmodel, matlitho, vars,verbosecpu); if (nfreq>1){ + EarthRheology(vars,frequencies_fluid,1,matlitho,femmodel); compute_love_numbers(Fluid, NULL, forcing_type, sh_nmax,frequencies_fluid, femmodel, matlitho, vars,verbosecpu); sh_cutoff=sh_nmax; @@ -1402,16 +1910,8 @@ } else sh_cutoff=sh_nmax; - - - - //Take care of rotationnal feedback love numbers first, if relevant /*{{{*/ - if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){ // if forcing is surface loading and we have degree 2 - if(VerboseSolution() && verbosecpu) _printf_(" tidal\n"); - int tidal_forcing_type=9; - compute_love_numbers(Tidalf_local, NULL,tidal_forcing_type=9, 2,frequencies_local, femmodel, matlitho, vars,verbosecpu); - } - /*}}}*/ - - //Resume requested forcing_type + + delete Fluid; + + //Requested forcing_type if (nfreq>1){ // if we are not running just elastic love numbers if(VerboseSolution() && verbosecpu){ @@ -1420,8 +1920,18 @@ else _printf_(" love\n"); } + EarthRheology(vars,frequencies_local,nf_local,matlitho,femmodel); compute_love_numbers(Lovef_local, Elastic, forcing_type, sh_cutoff, frequencies_local, femmodel, matlitho, vars,verbosecpu); } else{ Lovef_local->Copy(Elastic); + } + /*}}}*/ + + //Take care of rotationnal feedback love numbers, if relevant /*{{{*/ + if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){ // if forcing is surface loading and we have degree 2 + if(VerboseSolution() && verbosecpu) _printf_(" tidal\n"); + int tidal_forcing_type=9; + //no need to call EarthRheology, we already have the right one + compute_love_numbers(Tidalf_local, NULL,tidal_forcing_type=9, 2,frequencies_local, femmodel, matlitho, vars,verbosecpu); } /*}}}*/ @@ -1431,14 +1941,15 @@ if (istemporal && !complex_computation){ /*Initialize*/ - doubletype* pmtf_colineart=NULL; - doubletype* pmtf_orthot=NULL; - LoveNumbers* Lovet=NULL; - LoveNumbers* Tidalt=NULL; /*Downcast arrays to be exported in parameters*/ - LoveNumbers* LovetDouble=NULL; - LoveNumbers* TidaltDouble=NULL; IssmDouble* pmtf_colineartDouble=NULL; IssmDouble* pmtf_orthotDouble=NULL; + /* parallel computing */ + LoveNumbers* LovefDouble_local=NULL; + LoveNumbers* LovetDouble_local=NULL; + LoveNumbers* TidaltDouble_local=NULL; + IssmDouble* pmtf_colineartDouble_local=NULL; + IssmDouble* pmtf_orthotDouble_local=NULL; + doubletype* pmtf_colineart_local=NULL; doubletype* pmtf_orthot_local=NULL; @@ -1446,11 +1957,6 @@ LoveNumbers* Tidalt_local=NULL; - Lovet= new LoveNumbers(sh_nmin,sh_nmax,nt,matlitho); - Tidalt= new LoveNumbers(2,2,nt,matlitho); - pmtf_colineart=xNewZeroInit(3*nt); - pmtf_orthot=xNewZeroInit(3*nt); - - Lovet_local= new LoveNumbers(sh_nmin,sh_nmax,nt_local,matlitho); - Tidalt_local= new LoveNumbers(2,2,nt_local,matlitho); + Lovet_local= new LoveNumbers(sh_nmin,sh_nmax,nt_local,lower_row/2/NTit,nt,matlitho); + Tidalt_local= new LoveNumbers(2,2,nt_local,lower_row/2/NTit,nt,matlitho); pmtf_colineart_local=xNewZeroInit(3*nt_local); pmtf_orthot_local=xNewZeroInit(3*nt_local); @@ -1458,99 +1964,131 @@ love_freq_to_temporal(Lovet_local,Tidalt_local,pmtf_colineart_local,pmtf_orthot_local,Lovef_local,Tidalf_local,frequencies_local,femmodel,verbosecpu); - /* MPI Gather */ /*{{{*/ - Lovef->LoveMPI_Gather(Lovef_local, lower_row); - if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){ - Tidalf->LoveMPI_Gather(Tidalf_local, lower_row); - } - Lovet->LoveMPI_Gather(Lovet_local, lower_row/2/NTit); - Tidalt->LoveMPI_Gather(Tidalt_local, lower_row/2/NTit); + + if(VerboseSolution() && verbosecpu) _printf_(" Assembling parralel vectors..."); + + //delete Lovef_local; + delete Tidalf_local; + //Lovet + LovetDouble_local= new LoveNumbers(sh_nmin,sh_nmax,nt_local,lower_row/2/NTit,nt,matlitho); + Lovet_local->DownCastToDouble(LovetDouble_local); + delete Lovet_local; + LovetDouble_local->Broadcast(); + + //Lovef + LovefDouble_local= new LoveNumbers(sh_nmin,sh_nmax,nf_local,lower_row,nfreq,matlitho); + Lovef_local->DownCastToDouble(LovefDouble_local); + delete Lovef_local; + LovefDouble_local->Broadcast(); + + if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){ + TidaltDouble_local= new LoveNumbers(2,2,nt_local,lower_row,nfreq,matlitho); + Tidalt_local->DownCastToDouble(TidaltDouble_local); + delete Tidalt_local; + TidaltDouble_local->Broadcast(); + } + //pmtf: + pmtf_colineartDouble_local=xNew(nt_local); + pmtf_orthotDouble_local=xNew(nt_local); + /*Downcast*/ /*{{{*/ + for(int i=0;i(pmtf_colineart_local[i*3+2]); + pmtf_orthotDouble_local[i]=DownCastVarToDouble(pmtf_orthot_local[i*3+2]); + } + /*}}}*/ + xDelete(pmtf_colineart_local); + xDelete(pmtf_orthot_local); + pmtf_colineartDouble=xNew(nt); + pmtf_orthotDouble=xNew(nt); + int* recvcounts=xNew(IssmComm::GetSize()); int* displs=xNew(IssmComm::GetSize()); int rc; int offset; - rc=3*nt_local; - offset=3*lower_row/2/NTit; + rc=nt_local; + offset=lower_row/2/NTit; ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(pmtf_colineart_local, rc, ISSM_MPI_DOUBLE, pmtf_colineart, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(pmtf_orthot_local, rc, ISSM_MPI_DOUBLE, pmtf_orthot, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(pmtf_colineartDouble_local, rc, ISSM_MPI_DOUBLE, pmtf_colineartDouble, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(pmtf_orthotDouble_local, rc, ISSM_MPI_DOUBLE, pmtf_orthotDouble, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); xDelete(recvcounts); xDelete(displs); + + xDelete(pmtf_colineartDouble_local); + xDelete(pmtf_orthotDouble_local); /*}}}*/ - /*Downcast and add into parameters:*/ /*{{{*/ - LovetDouble= new LoveNumbers(sh_nmin,sh_nmax,nt,matlitho); - TidaltDouble= new LoveNumbers(2,2,nt,matlitho); - - pmtf_colineartDouble=xNew(nt); - pmtf_orthotDouble=xNew(nt); - - Lovet->DownCastToDouble(LovetDouble); - Tidalt->DownCastToDouble(TidaltDouble); - for(int i=0;iparameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LovetDouble->H,(sh_nmax+1)*nt,1)); - femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LovetDouble->K,(sh_nmax+1)*nt,1)); - femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LovetDouble->L,(sh_nmax+1)*nt,1)); - } + femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LovetDouble_local->H,(sh_nmax+1)*nt,1)); + femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LovetDouble_local->K,(sh_nmax+1)*nt,1)); + femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LovetDouble_local->L,(sh_nmax+1)*nt,1)); + } else if(forcing_type==11){ //surface loading - femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LovetDouble->H,(sh_nmax+1)*nt,1)); - femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LovetDouble->K,(sh_nmax+1)*nt,1)); - femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LovetDouble->L,(sh_nmax+1)*nt,1)); - femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,TidaltDouble->H,3*nt,1)); - femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,TidaltDouble->K,3*nt,1)); - femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,TidaltDouble->L,3*nt,1)); + femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LovetDouble_local->H,(sh_nmax+1)*nt,1)); + femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LovetDouble_local->K,(sh_nmax+1)*nt,1)); + femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LovetDouble_local->L,(sh_nmax+1)*nt,1)); + femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,TidaltDouble_local->H,3*nt,1)); + femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,TidaltDouble_local->K,3*nt,1)); + femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,TidaltDouble_local->L,3*nt,1)); femmodel->parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionColinearEnum,pmtf_colineartDouble,nt,1)); femmodel->parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionOrthogonalEnum,pmtf_orthotDouble,nt,1)); } + /*}}}*/ + + /*Add into external results*/ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKtEnum,LovetDouble_local->K,nt,sh_nmax+1,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHtEnum,LovetDouble_local->H,nt,sh_nmax+1,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLtEnum,LovetDouble_local->L,nt,sh_nmax+1,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKfEnum,LovefDouble_local->K,nfreq,sh_nmax+1,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHfEnum,LovefDouble_local->H,nfreq,sh_nmax+1,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLfEnum,LovefDouble_local->L,nfreq,sh_nmax+1,0,0)); + + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveTidalKtEnum,TidaltDouble_local->K,nt,3,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveTidalHtEnum,TidaltDouble_local->H,nt,3,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveTidalLtEnum,TidaltDouble_local->L,nt,3,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LovePMTF1tEnum,pmtf_colineartDouble,nt,1,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LovePMTF2tEnum,pmtf_orthotDouble,nt,1,0,0)); + /*Only when love_kernels is on*/ + if (love_kernels==1) { + // femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsEnum,LovefDouble_local->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0)); + } xDelete(pmtf_colineartDouble); xDelete(pmtf_orthotDouble); - /*}}}*/ - - /*Add into external results, no need to downcast, we can handle complexes/quad precision: */ - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKtEnum,Lovet->K,nt,sh_nmax+1,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHtEnum,Lovet->H,nt,sh_nmax+1,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLtEnum,Lovet->L,nt,sh_nmax+1,0,0)); - - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveTidalKtEnum,Tidalt->K,nt,3,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveTidalHtEnum,Tidalt->H,nt,3,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveTidalLtEnum,Tidalt->L,nt,3,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LovePMTF1tEnum,pmtf_colineart,nt,3,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LovePMTF2tEnum,pmtf_orthot,nt,3,0,0)); - - delete Lovet; - delete Tidalt; - delete Lovet_local; - delete Tidalt_local; - delete LovetDouble; - delete TidaltDouble; - - xDelete(pmtf_colineart); - xDelete(pmtf_orthot); + delete LovetDouble_local; + delete TidaltDouble_local; + } else{ LoveNumbers* LovefDouble=NULL; - LovefDouble= new LoveNumbers(sh_nmin,sh_nmax,nfreq,matlitho); + LoveNumbers* LovefDouble_local=NULL; + LovefDouble= new LoveNumbers(sh_nmin,sh_nmax,nfreq,lower_row,nfreq,matlitho); + LovefDouble_local= new LoveNumbers(sh_nmin,sh_nmax,nf_local,lower_row,nfreq,matlitho); + + LoveNumbers* TidalfDouble=NULL; + LoveNumbers* TidalfDouble_local=NULL; + TidalfDouble= new LoveNumbers(2,2,nfreq,lower_row,nfreq,matlitho); + TidalfDouble_local= new LoveNumbers(2,2,nf_local,lower_row,nfreq,matlitho); + + Lovef_local->DownCastToDouble(LovefDouble_local); + Tidalf_local->DownCastToDouble(TidalfDouble_local); /*MPI_Gather*/ if (nfreq>1){ - Lovef->LoveMPI_Gather(Lovef_local, lower_row); + LovefDouble->LoveMPI_Gather(LovefDouble_local, lower_row); if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){ - Tidalf->LoveMPI_Gather(Tidalf_local, lower_row); + TidalfDouble->LoveMPI_Gather(TidalfDouble_local, lower_row); } } else{ - Lovef->Copy(Elastic); - Tidalf->Copy(Tidalf_local); + Elastic->DownCastToDouble(LovefDouble); + Tidalf_local->DownCastToDouble(TidalfDouble); } /*Add into parameters:*/ - Lovef->DownCastToDouble(LovefDouble); if(forcing_type==9){ //tidal loading femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LovefDouble->H,(sh_nmax+1)*nfreq,1)); @@ -1563,14 +2101,33 @@ femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LovefDouble->L,(sh_nmax+1)*nfreq,1)); } + + /*Add into external results:*/ + if (complex_computation){ + //FIXME: complex external result not supported yet + //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKfEnum,Lovef->K,nfreq,sh_nmax+1,0,0)); + //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHfEnum,Lovef->H,nfreq,sh_nmax+1,0,0)); + //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLfEnum,Lovef->L,nfreq,sh_nmax+1,0,0)); + ///*Only when love_kernels is on*/ + //if (love_kernels==1) { + // femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsEnum,Lovef->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0)); + //} + } + else{ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKfEnum,LovefDouble->K,nfreq,sh_nmax+1,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHfEnum,LovefDouble->H,nfreq,sh_nmax+1,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLfEnum,LovefDouble->L,nfreq,sh_nmax+1,0,0)); + /*Only when love_kernels is on*/ + if (love_kernels==1) { + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsEnum,LovefDouble->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0)); + } + } + + + delete Lovef_local; delete LovefDouble; - } - - /*Add into external results:*/ - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKfEnum,Lovef->K,nfreq,sh_nmax+1,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHfEnum,Lovef->H,nfreq,sh_nmax+1,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLfEnum,Lovef->L,nfreq,sh_nmax+1,0,0)); - /*Only when love_kernels is on*/ - if (love_kernels==1) { - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsEnum,Lovef->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0)); + delete LovefDouble_local; + delete Tidalf_local; + delete TidalfDouble; + delete TidalfDouble_local; } /*Free resources:*/ @@ -1578,55 +2135,6 @@ xDelete(frequencies_local); xDelete(frequencies_elastic); - delete Lovef; - delete Lovef_local; - delete Tidalf; - delete Tidalf_local; + delete Elastic; - /* Legacy for fortran core, to be removed after complete validation */ - - //IssmDouble g0,r0; - //femmodel->parameters->FindParam(&g0,LoveG0Enum); - //femmodel->parameters->FindParam(&r0,LoveR0Enum); - //IssmDouble* LoveKr = xNewZeroInit(nfreq*(sh_nmax+1)); - //IssmDouble* LoveHr = xNewZeroInit(nfreq*(sh_nmax+1)); - //IssmDouble* LoveLr = xNewZeroInit(nfreq*(sh_nmax+1)); - //IssmDouble* LoveKi = xNewZeroInit(nfreq*(sh_nmax+1)); - //IssmDouble* LoveHi = xNewZeroInit(nfreq*(sh_nmax+1)); - //IssmDouble* LoveLi = xNewZeroInit(nfreq*(sh_nmax+1)); - - /*Initialize love kernels (real and imaginary parts): */ - //IssmDouble* LoveKernelsr= xNewZeroInit(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6); - //IssmDouble* LoveKernelsi= xNewZeroInit(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6); - - /*call the main module: */ - //if (false){ - //FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHr,LoveLr,LoveLi,LoveKernelsr,LoveKernelsi, //output - // nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs - // matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu, - // matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->rheologymodel, matlitho->issolid //matlitho inputs - // ); - - /*Only when love_kernels is on*/ - //if (love_kernels==1) { - //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernelsr,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0)); - //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0)); - //} - - /*Add love matrices to results:*/ - //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKtEnum,LoveKr,sh_nmax+1,nfreq,0,0)); - //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHtEnum,LoveHr,sh_nmax+1,nfreq,0,0)); - //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLtEnum,LoveLr,sh_nmax+1,nfreq,0,0)); - //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKfEnum,LoveKi,sh_nmax+1,nfreq,0,0)); - //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHfEnum,LoveHi,sh_nmax+1,nfreq,0,0)); - //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLfEnum,LoveLi,sh_nmax+1,nfreq,0,0)); - - //xDelete(LoveKr); - //xDelete(LoveHr); - //xDelete(LoveLr); - //xDelete(LoveKernelsr); - //xDelete(LoveKi); - //xDelete(LoveHi); - //xDelete(LoveLi); - //xDelete(LoveKernelsi); } /*}}}*/ @@ -1634,39 +2142,83 @@ /*cores and template instantiations:*/ /*template instantiations :{{{*/ +// IssmDouble template void love_core_template(FemModel* femmodel); -template void love_core_template(FemModel* femmodel); -template void fill_yi_prefactor(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); -template void fill_yi_prefactor(IssmComplex* yi_prefactor, int* pdeg, IssmComplex* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); -template void GetEarthRheology(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega, Matlitho* matlitho); -template void GetEarthRheology(IssmComplex* pla, IssmComplex* pmu, int layer_index, IssmComplex omega, Matlitho* matlitho); -template void yi_boundary_conditions(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type); -template void yi_boundary_conditions(IssmComplex* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type); -template void yi_derivatives(IssmDouble* dydx, IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); -template void yi_derivatives(IssmComplex* dydx, IssmComplex* y, int layer_index, int n, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); -template void propagate_yi_RK2(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); -template void propagate_yi_RK2(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); -template void propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); -template void propagate_yi_euler(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); -template void Innersphere_boundaryconditions(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); -template void Innersphere_boundaryconditions(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); -template void build_yi_system(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); -template void build_yi_system(IssmComplex* yi, int deg, IssmComplex omega, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); -template void solve_yi_system(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* frequencies, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars,bool verbosecpu); -template void solve_yi_system(IssmComplex* loveh, IssmComplex* lovel, IssmComplex* lovek, int deg, IssmComplex omega, IssmDouble* frequencies, IssmComplex* yi, IssmComplex* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars,bool verbosecpu); -template void compute_love_numbers(LoveNumbers* Lovef, LoveNumbers* Elastic, int forcing_type, int sh_cutoff,IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu); -template void compute_love_numbers(LoveNumbers* Lovef, LoveNumbers* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu); +template LoveVariables* love_init(FemModel* femmodel, Matlitho* matlitho,bool verbosecpu); +template void fill_yi_prefactor(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void GetEarthRheology(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega, Matlitho* matlitho, FemModel* femmodel); +template IssmDouble GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); +template void yi_boundary_conditions(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type); +template void yi_derivatives(IssmDouble* dydx, IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void propagate_yi_RK2(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void Innersphere_boundaryconditions(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void Coremantle_boundaryconditions(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void build_yi_system(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); +template void solve_yi_system(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* frequencies, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars,bool verbosecpu); +template void compute_love_numbers(LoveNumbers* Lovef, LoveNumbers* Elastic, int forcing_type, int sh_cutoff,IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu); template IssmDouble factorial(int n); template IssmDouble* postwidder_coef(int NTit); template IssmDouble n_C_r(int n, int r); template void postwidder_transform(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int sh_nmax,int NTit, IssmDouble* xi, FemModel* femmodel); +template void EarthRheology(LoveVariables* vars, IssmDouble* frequencies, int nfreq, Matlitho* matlitho, FemModel* femmodel); +template IssmDouble HypergeomTableLookup(IssmDouble z1, IssmDouble alpha, IssmDouble* h1, IssmDouble* z, int nz, int nalpha); + +//IssmComplex +template void love_core_template(FemModel* femmodel); +template LoveVariables* love_init(FemModel* femmodel, Matlitho* matlitho,bool verboscpu); +template void fill_yi_prefactor(IssmComplex* yi_prefactor, int* pdeg, IssmComplex* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void GetEarthRheology(IssmComplex* pla, IssmComplex* pmu, int layer_index, IssmComplex omega, Matlitho* matlitho, FemModel* femmodel); +template IssmComplex GetGravity(IssmComplex r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); +template void yi_boundary_conditions(IssmComplex* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type); +template void yi_derivatives(IssmComplex* dydx, IssmComplex* y, int layer_index, int n, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void propagate_yi_RK2(IssmComplex* y, IssmComplex xmin, IssmComplex xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void propagate_yi_RK4(IssmComplex* y, IssmComplex xmin, IssmComplex xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void propagate_yi_euler(IssmComplex* y, IssmComplex xmin, IssmComplex xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void Innersphere_boundaryconditions(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void Coremantle_boundaryconditions(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void build_yi_system(IssmComplex* yi, int deg, IssmComplex omega, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); +template void solve_yi_system(IssmComplex* loveh, IssmComplex* lovel, IssmComplex* lovek, int deg, IssmComplex omega, IssmDouble* frequencies, IssmComplex* yi, IssmComplex* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars,bool verbosecpu); +template void compute_love_numbers(LoveNumbers* Lovef, LoveNumbers* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu); +template void EarthRheology(LoveVariables* vars, IssmDouble* frequencies, int nfreq, Matlitho* matlitho, FemModel* femmodel); +template IssmComplex HypergeomTableLookup(IssmComplex z1, IssmComplex alpha, IssmDouble* h1, IssmDouble* z, int nz, int nalpha); + +//__float128 +#ifdef _HAVE_MPLAPACK_ +template void love_core_template<__float128>(FemModel* femmodel); +template LoveVariables<__float128>* love_init<__float128>(FemModel* femmodel, Matlitho* matlitho, bool verbosecpu); +template void fill_yi_prefactor<__float128>(__float128* yi_prefactor, int* pdeg, __float128* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars); +template void GetEarthRheology<__float128>(__float128* pla, __float128* pmu, int layer_index, __float128 omega, Matlitho* matlitho, FemModel* femmodel); +template __float128 GetGravity<__float128>(__float128 r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables<__float128>* vars); +template void yi_boundary_conditions<__float128>(__float128* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables<__float128>* vars, int forcing_type); +template void yi_derivatives<__float128>(__float128* dydx, __float128* y, int layer_index, int n, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars); +template void propagate_yi_RK2<__float128>(__float128* y, __float128 xmin, __float128 xmax, int layer_index, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars); +template void propagate_yi_RK4<__float128>(__float128* y, __float128 xmin, __float128 xmax, int layer_index, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars); +template void propagate_yi_euler<__float128>(__float128* y, __float128 xmin, __float128 xmax, int layer_index, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars); +template void Innersphere_boundaryconditions<__float128>(__float128* yi, int layer_index, int deg, __float128 omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars); +template void Coremantle_boundaryconditions<__float128>(__float128* yi, int layer_index, int deg, __float128 omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars); +template void build_yi_system<__float128>(__float128* yi, int deg, __float128 omega, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables<__float128>* vars); +template void solve_yi_system<__float128>(__float128* loveh, __float128* lovel, __float128* lovek, int deg, __float128 omega, IssmDouble* frequencies, __float128* yi, __float128* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars,bool verbosecpu); +template void compute_love_numbers<__float128>(LoveNumbers<__float128>* Lovef, LoveNumbers<__float128>* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars, bool verbosecpu); +template __float128 factorial<__float128>(int n); +template __float128* postwidder_coef<__float128>(int NTit); +template __float128 n_C_r<__float128>(int n, int r); +template void postwidder_transform<__float128>(__float128* Lovet, __float128* Lovef,int d, int t, int sh_nmax,int NTit, __float128* xi, FemModel* femmodel); +template void EarthRheology<__float128>(LoveVariables<__float128>* vars, IssmDouble* frequencies, int nfreq, Matlitho* matlitho, FemModel* femmodel); +template __float128 HypergeomTableLookup(__float128 z1, __float128 alpha, IssmDouble* h1, IssmDouble* z, int nz, int nalpha); +#endif /*}}}*/ void love_core(FemModel* femmodel){ /*{{{*/ - bool complex_computation; + bool quad_precision; femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum); + femmodel->parameters->FindParam(&quad_precision,LoveQuadPrecisionEnum); if(complex_computation) love_core_template(femmodel); +#ifdef _HAVE_MPLAPACK_ + else if(quad_precision) love_core_template<__float128>(femmodel); +#endif else love_core_template(femmodel); Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 27308) @@ -22,5 +22,5 @@ void TransferForcing(FemModel* femmodel,int forcingenum); void TransferSealevel(FemModel* femmodel,int forcingenum); -bool slcconvergence(Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs); +bool slcconvergence(IssmDouble* RSLg,IssmDouble* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs, IssmDouble totaloceanarea,FemModel* femmodel); IssmDouble SealevelloadsOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble totaloceanarea); void PolarMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom, bool computefuture); @@ -52,4 +52,7 @@ femmodel->SetCurrentConfiguration(SealevelchangeAnalysisEnum); + /*Run coupler input transfer:*/ + couplerinput_core(femmodel); + /*run geometry core: */ slgeom=sealevelchange_geometry(femmodel); @@ -57,7 +60,4 @@ /*any external forcings?:*/ solidearthexternal_core(femmodel); - - /*Run coupler input transfer:*/ - couplerinput_core(femmodel); /*Run geodetic:*/ @@ -184,9 +184,12 @@ /*if we are carrying loads but are not yet computing grd core, accumulate them and skip * the rest: */ - if (countparameters->SetParam(count,SealevelchangeRunCountEnum); + if (count!=frequency){ + if (count>frequency){ + count=1; + femmodel->parameters->SetParam(count,SealevelchangeRunCountEnum); + } return; } + /*Basins are supposed to accumulate loads and hand them over to the Earth @@ -208,4 +211,6 @@ TransferForcing(femmodel,DeltaBottomPressureEnum); TransferForcing(femmodel,DeltaTwsEnum); + TransferForcing(femmodel,MaskOceanLevelsetEnum); + TransferForcing(femmodel,MaskIceLevelsetEnum); /*transfer external forcings back to Earth:*/ @@ -228,5 +233,5 @@ GrdLoads* loads=NULL; - Vector* oldsealevelloads=NULL; + IssmDouble* oldsealevelloads=NULL; Vector* oceanareas=NULL; IssmDouble totaloceanarea; @@ -249,4 +254,5 @@ int grdmodel; int sealevelloading=0; + bool sal=false; bool viscous=false; bool rotation=false; @@ -256,7 +262,4 @@ /*}}}*/ - - /*Verbose: */ - if(VerboseSolution()) _printf0_(" computing GRD patterns\n"); /*retrieve parameters:{{{*/ @@ -290,4 +293,8 @@ } + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing GRD patterns\n"); + + /*retrieve parameters: {{{*/ femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); @@ -298,4 +305,5 @@ femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum); femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + femmodel->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum); /*}}}*/ @@ -336,7 +344,12 @@ PolarMotion(&polarmotionvector[0],femmodel,loads,slgeom,computefuture=false); - oldsealevelloads=loads->vsealevelloads->Duplicate(); loads->vsealevelloads->Copy(oldsealevelloads); + oldsealevelloads=xNewZeroInit(nel); + if (loads->sealevelloads){ + xMemCpy(oldsealevelloads,loads->sealevelloads,nel); + } + /*convolve load and sealevel loads on oceans:*/ + loads->Combineloads(nel,slgeom); //This combines loads and sealevelloads into a single vector for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); @@ -367,15 +380,17 @@ loads->BroadcastSealevelLoads(); + if (!sal) xDelete(oldsealevelloads); break; + //convergence? - if(slcconvergence(loads->vsealevelloads,oldsealevelloads,eps_rel,eps_abs)){ - delete oldsealevelloads; break; + if(slcconvergence(loads->sealevelloads,oldsealevelloads,eps_rel,eps_abs,totaloceanarea,femmodel)){ + xDelete(oldsealevelloads); break; } //early return? if(iterations>=max_nonlinear_iterations){ - delete oldsealevelloads; break; + xDelete(oldsealevelloads); break; } iterations++; - delete oldsealevelloads; + xDelete(oldsealevelloads); } @@ -388,4 +403,5 @@ /*convolve loads and sea level loads to get the deformation:*/ + loads->Combineloads(nel,slgeom); //This combines loads and sealevelloads into a single vector for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); @@ -512,9 +528,15 @@ int iscoupling; int horiz=0; + int count, frequency; /*retrieve more parameters:*/ femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum); femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - + femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); + femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); + + count++; + femmodel->parameters->SetParam(count,SealevelchangeRunCountEnum); + if(iscoupling){ /*transfer sea level back to ice caps:*/ @@ -526,4 +548,5 @@ } } + }; /*}}}*/ void ivins_deformation_core(FemModel* femmodel){ /*{{{*/ @@ -588,5 +611,6 @@ IssmDouble* areae = NULL; int nel; - int* lids; + int* lids=NULL; + int* n_activevertices=NULL; int grdmodel=0; @@ -596,5 +620,5 @@ /*early return?:*/ - if(grdmodel==IvinsEnum) return; + if(grdmodel!=ElasticEnum) return; /*Verbose: */ @@ -607,8 +631,17 @@ /*Compute element ids, used to speed up computations in convolution phase:{{{*/ lids=xNew(femmodel->vertices->Size()); + n_activevertices = xNew(nel); + //initialize lids to -1, vertex count to 3 + for (int v=0; vvertices->Size();v++) lids[v]=-1; + for (int e=0; eelements->objects){ Element* element=xDynamicCast(object); for(int i=0;i<3;i++){ + // if lids where we are looking points to an element id (.i.e. not -1) then we are about to claim that element's vertex + // and need to lower the number of vertices it is in charge of + if (lids[element->vertices[i]->lid] !=-1){ + n_activevertices[lids[element->vertices[i]->lid]]-=1; + } lids[element->vertices[i]->lid]=element->lid; } @@ -620,5 +653,5 @@ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - element->SealevelchangeGeometryInitial(xxe,yye,zze,areae,lids); + element->SealevelchangeGeometryInitial(xxe,yye,zze,areae,lids,n_activevertices); } @@ -639,5 +672,6 @@ xDelete(zze); xDelete(areae); - xDelete(lids); + xDelete(lids); + xDelete(n_activevertices); return; @@ -658,4 +692,5 @@ int grdmodel=0; int isgrd=0; + int count, frequency; SealevelGeometry* slgeom=NULL; @@ -663,5 +698,8 @@ femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum); + femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); + femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); if(grdmodel!=ElasticEnum || !isgrd) return NULL; + if(count!=frequency)return NULL; /*retrieve parameters:*/ @@ -722,4 +760,5 @@ int isgrd=0; int horiz=0; + int count, frequency; /*early return?:*/ @@ -727,5 +766,8 @@ femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum); femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); + femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); if(grdmodel!=ElasticEnum || !isgrd) return; + if(count!=frequency)return; for (int l=0;l* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ - +bool slcconvergence(IssmDouble* RSLg,IssmDouble* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs, IssmDouble totaloceanarea, FemModel* femmodel){ /*{{{*/ + + int nel; bool converged=true; - IssmDouble ndS,nS; - Vector *dRSLg = NULL; + IssmDouble ndS,nS, nS_old; + IssmDouble* dRSLg = NULL; + IssmDouble rho_water =0; + + femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum); + femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); //compute norm(du) and norm(u) if requested - dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0); - ndS=dRSLg->Norm(NORM_TWO); + dRSLg=xNewZeroInit(nel); + + ndS=0; + nS=0; + nS_old=0; + + for (int e=0;e(ndS)){ - _error_("convergence criterion is NaN (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")"); + _error_("convergence criterion is NaN (RSL_old=" << nS_old << " RSL=" << nS << ")"); } if(!xIsNan(eps_rel)){ - nS=RSLg_old->Norm(NORM_TWO); - if (xIsNan(nS)) _error_("convergence criterion is NaN! (check the initial RSL)"); + if (xIsNan(nS_old)) _error_("convergence criterion is NaN! (check the initial RSL)"); } //clean up - delete dRSLg; + xDelete(dRSLg); //print @@ -1048,4 +1105,5 @@ IssmDouble* forcing=NULL; Vector* forcingglobal=NULL; + IssmDouble* transfercount=NULL; int* nvs=NULL; @@ -1088,5 +1146,9 @@ nv=femmodel->vertices->NumberOfVertices(); existforcing=reCast(femmodel->inputs->Exist(forcingenum)); - if(existforcing)GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum); + if(existforcing){ + GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum); + GetVectorFromInputsx(&transfercount,femmodel,CouplingTransferCountEnum,VertexSIdEnum); + for (int i=0;iSet(0.0); + /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/ femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelchangeTransitionsEnum); @@ -1146,4 +1210,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,ADD_VAL); xDelete(index); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27307) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27308) @@ -312,4 +312,5 @@ LoveMinIntegrationStepsEnum, LoveMaxIntegrationdrEnum, + LoveIntegrationSchemeEnum, LoveKernelsEnum, LoveMu0Enum, @@ -935,4 +936,5 @@ SealevelUNorthEsaEnum, SealevelchangeIndicesEnum, + SealevelchangeConvolutionVerticesEnum, SealevelchangeAlphaIndexEnum, SealevelchangeAzimuthIndexEnum, @@ -953,4 +955,5 @@ SealevelchangeViscousNEnum, SealevelchangeViscousEEnum, + CouplingTransferCountEnum, SedimentHeadEnum, SedimentHeadOldEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27308) @@ -320,4 +320,5 @@ case LoveMinIntegrationStepsEnum : return "LoveMinIntegrationSteps"; case LoveMaxIntegrationdrEnum : return "LoveMaxIntegrationdr"; + case LoveIntegrationSchemeEnum : return "LoveIntegrationScheme"; case LoveKernelsEnum : return "LoveKernels"; case LoveMu0Enum : return "LoveMu0"; @@ -941,4 +942,5 @@ case SealevelUNorthEsaEnum : return "SealevelUNorthEsa"; case SealevelchangeIndicesEnum : return "SealevelchangeIndices"; + case SealevelchangeConvolutionVerticesEnum : return "SealevelchangeConvolutionVertices"; case SealevelchangeAlphaIndexEnum : return "SealevelchangeAlphaIndex"; case SealevelchangeAzimuthIndexEnum : return "SealevelchangeAzimuthIndex"; @@ -959,4 +961,5 @@ case SealevelchangeViscousNEnum : return "SealevelchangeViscousN"; case SealevelchangeViscousEEnum : return "SealevelchangeViscousE"; + case CouplingTransferCountEnum : return "CouplingTransferCount"; case SedimentHeadEnum : return "SedimentHead"; case SedimentHeadOldEnum : return "SedimentHeadOld"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27307) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27308) @@ -326,4 +326,5 @@ else if (strcmp(name,"LoveMinIntegrationSteps")==0) return LoveMinIntegrationStepsEnum; else if (strcmp(name,"LoveMaxIntegrationdr")==0) return LoveMaxIntegrationdrEnum; + else if (strcmp(name,"LoveIntegrationScheme")==0) return LoveIntegrationSchemeEnum; else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; @@ -382,9 +383,9 @@ else if (strcmp(name,"OutputBufferSizePointer")==0) return OutputBufferSizePointerEnum; else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; - else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; + if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; + else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; else if (strcmp(name,"SmbARMAInitialTime")==0) return SmbARMAInitialTimeEnum; - else if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum; + if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum; + else if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum; else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum; else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; - else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; + if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; + else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; - else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; + if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; + else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum; else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum; - else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum; + if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum; + else if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum; else if (strcmp(name,"MaterialsRheologyEbar")==0) return MaterialsRheologyEbarEnum; else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum; @@ -962,4 +963,5 @@ else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum; else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum; + else if (strcmp(name,"SealevelchangeConvolutionVertices")==0) return SealevelchangeConvolutionVerticesEnum; else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum; else if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum; @@ -980,4 +982,5 @@ else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum; + else if (strcmp(name,"CouplingTransferCount")==0) return CouplingTransferCountEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; @@ -995,11 +998,11 @@ else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum; else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum; - else if (strcmp(name,"SmbA")==0) return SmbAEnum; - else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum; - else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum; + if (strcmp(name,"SmbA")==0) return SmbAEnum; + else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum; + else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum; + else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum; else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum; else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; @@ -1118,11 +1121,11 @@ else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; - else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum; - else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; - else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; + if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum; + else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; + else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; + else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; @@ -1241,11 +1244,11 @@ else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; - else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum; - else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; - else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; + if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum; + else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; + else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum; + else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum; else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum; @@ -1364,11 +1367,11 @@ else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; - else if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum; - else if (strcmp(name,"DebrisSolution")==0) return DebrisSolutionEnum; - else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; + if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum; + else if (strcmp(name,"DebrisSolution")==0) return DebrisSolutionEnum; + else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; + else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; else if (strcmp(name,"Dense")==0) return DenseEnum; else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum; @@ -1487,11 +1490,11 @@ else if (strcmp(name,"LoveHf")==0) return LoveHfEnum; else if (strcmp(name,"LoveHt")==0) return LoveHtEnum; - else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; - else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; - else if (strcmp(name,"LoveKf")==0) return LoveKfEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"LoveKt")==0) return LoveKtEnum; + if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; + else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; + else if (strcmp(name,"LoveKf")==0) return LoveKfEnum; + else if (strcmp(name,"LoveKt")==0) return LoveKtEnum; else if (strcmp(name,"LoveLf")==0) return LoveLfEnum; else if (strcmp(name,"LoveLt")==0) return LoveLtEnum; @@ -1610,11 +1613,11 @@ else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; - else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; - else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; - else if (strcmp(name,"Scaled")==0) return ScaledEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; + if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; + else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; + else if (strcmp(name,"Scaled")==0) return ScaledEnum; + else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum; else if (strcmp(name,"SealevelchangePolarMotionX")==0) return SealevelchangePolarMotionXEnum; Index: /issm/trunk-jpl/src/m/classes/love.m =================================================================== --- /issm/trunk-jpl/src/m/classes/love.m (revision 27308) +++ /issm/trunk-jpl/src/m/classes/love.m (revision 27308) @@ -0,0 +1,237 @@ +%love class definition +% +% Usage: +% md.love=love(); + +classdef love + properties (SetAccess=public) + nfreq = 0; + frequencies = 0; + sh_nmax = 0; + sh_nmin = 0; + g0 = 0; + r0 = 0; + mu0 = 0; + Gravitational_Constant = 0; + chandler_wobble = 0; + allow_layer_deletion = 0; + underflow_tol = 0; + pw_threshold = 0; + min_integration_steps = 0; + max_integration_dr = 0; + integration_scheme = 0; + istemporal = 0; + n_temporal_iterations = 0; + time = 0; + love_kernels = 0; + forcing_type = 0; + inner_core_boundary = 0; + core_mantle_boundary = 0; + complex_computation = 0; + quad_precision = 0; + debug=0; + hypergeom_table1=0; + hypergeom_table2=0; + hypergeom_nalpha=0; + hypergeom_nz=0; + hypergeom_z=0; + end + methods (Static) + function self = loadobj(self) % {{{ + % This function is directly called by matlab when a model object is + % loaded. Update old properties here + self.min_integration_steps = self.int_steps_per_layers; + end% }}} + end + methods + function self = love(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + %we setup an elastic love number computation by default. + self.nfreq=1; + self.frequencies=[0]; %Hz + self.sh_nmax=256; % .35 degree, 40 km at the equator. + self.sh_nmin=1; + % work on matlab script for computing g0 for given Earth's structure. + self.g0=9.81; % m/s^2; + self.r0=6371*1e3; %m; + self.mu0=1e11; % Pa + self.Gravitational_Constant=6.67259e-11; % m^3 kg^-1 s^-2 + self.chandler_wobble=0; + self.allow_layer_deletion=1; + self.underflow_tol=1e-16; %threshold of deep to surface love number ratio to trigger the deletion of layer + self.pw_threshold=1e-3; %if relative variation across frequencies is smaller than this ratio, the post-widder transform for time-dependent love numbers is bypassed + self.min_integration_steps = 50; + self.max_integration_dr = 10e3; + self.integration_scheme=1; + self.istemporal=0; + self.n_temporal_iterations=8; + self.time=[0]; %s + self.love_kernels=0; + self.forcing_type = 11; % surface loading + self.inner_core_boundary=1; + self.core_mantle_boundary=2; + self.complex_computation=0; + self.quad_precision=0; + self.debug=0; + self.hypergeom_table1=1; + self.hypergeom_table2=1; + self.hypergeom_nalpha=1; + self.hypergeom_nz=1; + self.hypergeom_z=0; + end % }}} + function disp(self) % {{{ + fielddisplay(self,'nfreq','number of frequencies sampled (default: 1, elastic) [Hz]'); + fielddisplay(self,'frequencies','frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'); + fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)'); + fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default: 1)'); + fielddisplay(self,'g0','adimensioning constant for gravity (default: 10) [m/s^2]'); + fielddisplay(self,'r0','adimensioning constant for radius (default: 6371*10^3) [m]'); + fielddisplay(self,'mu0','adimensioning constant for stress (default: 10^11) [Pa]'); + fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])'); + fielddisplay(self,'chandler_wobble','includes the inertial terms for the chandler wobble in the rotational feedback love numbers, only for forcing_type=11 (default: 0) (/!\ 1 is untested yet)'); + fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'); + fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)'); + fielddisplay(self,'pw_threshold','if relative variation across frequencies is smaller than this ratio, the post-widder transform for time-dependent love numbers is bypassed (default (1e-3)'); + fielddisplay(self,'min_integration_steps','minimum number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 50)'); + fielddisplay(self,'max_integration_dr','maximum length of radial steps to propagate the yi system from the bottom to the top of each layer (default: 10e3) [m]'); + fielddisplay(self,'integration_scheme','0: Euler, 1: Midpoint aka Runge-Kutta 2, 2: Runge-Kutta 4 (default: 1)'); + fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use love class function build_frequencies_from_time to meet consistency'}); + fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)'); + fielddisplay(self,'time','time vector for deformation if istemporal (default: 0) [s]'); + fielddisplay(self,'love_kernels','compute love numbers at depth? (default: 0)'); + fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):','1: Inner core boundary -- Volumic Potential','2: Inner core boundary -- Pressure','3: Inner core boundary -- Loading','4: Inner core boundary -- Tangential traction','5: Core mantle boundary -- Volumic Potential','6: Core mantle boundary -- Pressure','7: Core mantle boundary -- Loading','8: Core mantle boundary -- Tangential traction','9: Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '}); + fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)'); + fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)'); + fielddisplay(self,'complex_computation','return love numbers as 0: real (useful for elastic or temporal forms), 1: complex numbers (useful for Fourier spectral form) (default: 0)'); + fielddisplay(self,'quad_precision','toogle computation love numbers and post-widder transform with 32 digit precision, useful for temporal form (default: 1)'); + fielddisplay(self,'debug','outputs yi system matrix prior to solving (default: 0)'); + fielddisplay(self,'hypergeom_table1','table 1 for hypergeometric function, only for EBM rheology (default: [1])'); + fielddisplay(self,'hypergeom_table2','table 2 for hypergeometric function, only for EBM rheology (default: [1])'); + fielddisplay(self,'hypergeom_nalpha','length of hypergeometric table, only for EBM rheology (default: 1)'); + fielddisplay(self,'hypergeom_nz','width of hypergeometric table, only for EBM rheology (default: 1)'); + fielddisplay(self,'hypergeom_z','abscissa for hypergeometric table, only for EBM rheology (default: [0])'); + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ~ismember('LoveAnalysis',analyses), return; end + + md = checkfield(md,'fieldname','love.nfreq','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.frequencies','NaN',1,'Inf',1,'numel',md.love.nfreq); + md = checkfield(md,'fieldname','love.sh_nmax','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.sh_nmin','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.g0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.r0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.mu0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.Gravitational_Constant','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.chandler_wobble','values',[0 1]); + md = checkfield(md,'fieldname','love.allow_layer_deletion','values',[0 1]); + md = checkfield(md,'fieldname','love.underflow_tol','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.pw_threshold','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.min_integration_steps','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.max_integration_dr','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.integration_scheme','NaN',1,'Inf',1,'numel',1,'>=',0,'<=',2); + md = checkfield(md,'fieldname','love.love_kernels','values',[0 1]); + md = checkfield(md,'fieldname','love.forcing_type','NaN',1,'Inf',1,'numel',1,'>',0, '<=', 12); + md = checkfield(md,'fieldname','love.complex_computation','NaN',1,'Inf',1,'numel',1,'values',[0 1]); + md = checkfield(md,'fieldname','love.quad_precision','NaN',1,'Inf',1,'numel',1,'values',[0 1]); + md = checkfield(md,'fieldname','love.debug','NaN',1,'Inf',1,'numel',1,'values',[0 1]); + + md = checkfield(md,'fieldname','love.istemporal','values',[0 1]); + + md = checkfield(md,'fieldname','love.hypergeom_nalpha','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.hypergeom_nz','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.hypergeom_z','NaN',1,'Inf',1,'numel',md.love.hypergeom_nz); + md = checkfield(md,'fieldname','love.hypergeom_table1','NaN',1,'Inf',1,'numel',md.love.hypergeom_nz*md.love.hypergeom_nalpha); + md = checkfield(md,'fieldname','love.hypergeom_table2','NaN',1,'Inf',1,'numel',md.love.hypergeom_nz*md.love.hypergeom_nalpha); + + if md.love.istemporal==1 + md = checkfield(md,'fieldname','love.n_temporal_iterations','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.time','NaN',1,'Inf',1,'numel',md.love.nfreq/2/md.love.n_temporal_iterations); + end + if md.love.sh_nmin<=1 & (md.love.forcing_type==1 || md.love.forcing_type==5 || md.love.forcing_type==9) + error(['Degree 1 not supported for forcing type ' num2str(md.love.forcing_type) '. Use sh_min>=2 for this kind of calculation.']) + end + + if md.love.chandler_wobble==1 + disp('Warning, Chandler Wobble in Love number calculator has not been validated yet'); + end + + %need 'litho' material: + if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho')) + error('Need a ''litho'' material to run a Fourier Love number analysis'); + end + + mat=find(strcmpi(md.materials.nature,'litho')); + if md.love.forcing_type<=4 + md = checkfield(md,'fieldname','love.inner_core_boundary','NaN',1,'Inf',1,'numel',1,'>',0, '<=', md.materials(mat).numlayers); + elseif md.love.forcing_type<=8 + md = checkfield(md,'fieldname','love.core_mantle_boundary','NaN',1,'Inf',1,'numel',1,'>',0, '<=', md.materials(mat).numlayers); + end + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + WriteData(fid,prefix,'object',self,'fieldname','nfreq','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','frequencies','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'fieldname','sh_nmax','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','sh_nmin','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','g0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','r0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','mu0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','Gravitational_Constant','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','chandler_wobble','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','allow_layer_deletion','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','underflow_tol','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','pw_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','min_integration_steps','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','max_integration_dr','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','integration_scheme','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','istemporal','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','n_temporal_iterations','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','complex_computation','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','quad_precision','format','Boolean'); + %note: no need to marshall the time vector, we have frequencies + WriteData(fid,prefix,'object',self,'fieldname','love_kernels','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','forcing_type','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','inner_core_boundary','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','core_mantle_boundary','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','debug','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','hypergeom_table1','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','hypergeom_table2','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','hypergeom_nalpha','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','hypergeom_nz','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','hypergeom_z','format','DoubleMat','mattype',3); + + end % }}} + function self = extrude(self,md) % {{{ + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + error('not implemented yet!'); + end % }}} + function self=build_frequencies_from_time(self) % {{{ + if ~self.istemporal + error('cannot build frequencies for temporal love numbers if love.istemporal==0') + end + disp('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies'); + self.nfreq=length(self.time)*2*self.n_temporal_iterations; + self.frequencies=zeros(self.nfreq,1); + for i=1:length(self.time) + for j=1:2*self.n_temporal_iterations + if self.time(i)==0 + self.frequencies((i-1)*2*self.n_temporal_iterations +j) =0; % convention to avoid marshalling infinite numbers + else + self.frequencies((i-1)*2*self.n_temporal_iterations +j) = j*log(2)/self.time(i)/2/pi; + end + end + end + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 27307) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 27308) @@ -137,5 +137,5 @@ if isa(md.amr,'double'); md.amr=amr(); end %2017 Aug 29th - if isa(md.love,'double'); md.love=fourierlove(); end + if isa(md.love,'double'); md.love=love(); end %2017 Oct 26th if isa(md.calving,'calvingdev') @@ -289,5 +289,5 @@ md.calving = calving(); md.frontalforcings = frontalforcings(); - md.love = fourierlove(); + md.love = love(); md.esa = esa(); md.sampling = sampling(); Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 27307) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 27308) @@ -48,12 +48,18 @@ %is the coupler turned on? - for i=1:length(slm.icecaps), - if slm.icecaps{i}.transient.iscoupler==0, - warning(sprintf('sealevelmodel.m::checkconsistency: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name)); - end - end + %for i=1:length(slm.icecaps), + % if slm.icecaps{i}.transient.iscoupler==0, + % warning(sprintf('sealevelmodel.m::checkconsistency: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name)); + % end + %end - if slm.earth.transient.iscoupler==0, - warning('sealevelmodel.m::checkconsistency: earth model should have the transient coupler option turned on!'); + %if slm.earth.transient.iscoupler==0, + % warning('sealevelmodel.m::checkconsistency: earth model should have the transient coupler option turned on!'); + %end + + %check that the transition vectors have the right size: + + if slm.earth.mesh.numberofvertices ~= length(slm.earth.solidearth.transfercount) + error('sealevelmodel.m::checkconsistency: earth.solidearth.transfercount should be of size earth.mesh.numberofvertices') end @@ -75,5 +81,5 @@ for i=1:length(slm.icecaps), md= slm.icecaps{i}; - if ~isempty(find(md.dsl.steric_rate - slm.earth.dsl.steric_rate(slm.earth.dsl.transitions{i}))), + if ~isempty(find(md.dsl.sea_surface_height_above_geoid - slm.earth.dsl.sea_surface_height_above_geoid(slm.transitions{i}))), error(sprintf('sealevelmodel.m::checkconsistency: steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name)); end @@ -83,5 +89,5 @@ for i=1:length(slm.icecaps), md= slm.icecaps{i}; - if md.solidearthsettings.isgrd~=slm.earth.solidearthsettings.isgrd + if md.solidearth.settings.isgrd~=slm.earth.solidearth.settings.isgrd error(sprintf('sealevelmodel.m::checkconsistency: isgrd on ice cap %s is not the same as for the earth\n',md.miscellaneous.name)); end @@ -228,4 +234,5 @@ self.transitions={}; self.eltransitions={}; + self.earth.solidearth.transfercount=zeros(self.earth.mesh.numberofvertices,1); %for elements: @@ -248,4 +255,10 @@ self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force); + + self.earth.solidearth.transfercount(self.transitions{i})=self.earth.solidearth.transfercount(self.transitions{i})+1; + end + + for i=1:length(self.icecaps), + self.icecaps{i}.solidearth.transfercount=self.earth.solidearth.transfercount(self.transitions{i}); end end % }}} @@ -255,5 +268,5 @@ flags(self.transitions{i})=i; end - plotmodel(self.earth,'data',flags,'coastline','on'); + plotmodel(self.earth,'data',flags,'coastlines','on'); end % }}} Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 27307) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 27308) @@ -14,4 +14,5 @@ requested_outputs = {}; transitions = {}; + transfercount = []; partitionice = []; partitionhydro = []; @@ -51,4 +52,5 @@ fielddisplay(self,'planetradius','planet radius [m]'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); + fielddisplay(self,'transfercount','number of icecaps vertices are part of'); fielddisplay(self,'requested_outputs','additional outputs requested'); fielddisplay(self,'partitionice','ice partition vector for barystatic contribution'); @@ -71,4 +73,5 @@ %transitions should be a cell array of vectors: self.transitions={}; + self.transfercount=[0]; %no partitions requested for barystatic contribution: @@ -110,4 +113,7 @@ WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray'); + %WriteData(fid,prefix,'object',self,'fieldname','transfercount','data',self.transfercount,'format','DoubleMat','name', 'md.solidearth.transfercount','mattype',1); + %WriteData(fid,prefix,'data',self.transfercount,'format','DoubleMat','name', 'md.solidearth.transfercount'); + WriteData(fid,prefix,'object',self,'fieldname','transfercount','format','DoubleMat','mattype',1); if ~isempty(self.partitionice), Index: /issm/trunk-jpl/src/m/solve/solveslm.m =================================================================== --- /issm/trunk-jpl/src/m/solve/solveslm.m (revision 27307) +++ /issm/trunk-jpl/src/m/solve/solveslm.m (revision 27308) @@ -68,5 +68,5 @@ nps{end+1}=slm.earth.cluster.np; -BuildQueueScriptMultipleModels(cluster,slm.private.runtimename,slm.miscellaneous.name,slm.private.solution,valgrind,privateruntimenames,miscellaneousnames,nps); +BuildQueueScriptMultipleModels(cluster,slm.private.runtimename,slm.miscellaneous.name,slm.private.solution,privateruntimenames,miscellaneousnames,nps); %Upload all required files, given that each individual solution for icecaps and earth model already did: Index: /issm/trunk-jpl/test/NightlyRun/test2013.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2013.m (revision 27307) +++ /issm/trunk-jpl/test/NightlyRun/test2013.m (revision 27308) @@ -1,5 +1,5 @@ %Test Name: EarthSlc_Geometry -step=[1 2]; +step=[1]; if any(step==1) %mesh earth: Index: /issm/trunk-jpl/test/NightlyRun/test2070.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2070.m (revision 27307) +++ /issm/trunk-jpl/test/NightlyRun/test2070.m (revision 27308) @@ -40,5 +40,5 @@ md.love.pw_threshold=1e-3; md.love.Gravitational_Constant=6.6732e-11; -md.love.integration_steps_per_layer=100; +md.love.min_integration_steps=100; md.love.allow_layer_deletion=1; md.love.forcing_type=11; Index: /issm/trunk-jpl/test/NightlyRun/test2071.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2071.m (revision 27307) +++ /issm/trunk-jpl/test/NightlyRun/test2071.m (revision 27308) @@ -40,5 +40,5 @@ md.love.pw_threshold=1e-3; md.love.Gravitational_Constant=6.6732e-11; -md.love.integration_steps_per_layer=100; +md.love.min_integration_steps=100; md.love.allow_layer_deletion=1; md.love.forcing_type=11; Index: /issm/trunk-jpl/test/NightlyRun/test2072.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2072.m (revision 27307) +++ /issm/trunk-jpl/test/NightlyRun/test2072.m (revision 27308) @@ -44,5 +44,5 @@ md.love.pw_threshold=1e-3; md.love.Gravitational_Constant=6.6732e-11; -md.love.integration_steps_per_layer=100; +md.love.min_integration_steps=100; md.love.allow_layer_deletion=1; md.love.forcing_type=11; Index: /issm/trunk-jpl/test/NightlyRun/test2073.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2073.m (revision 27308) +++ /issm/trunk-jpl/test/NightlyRun/test2073.m (revision 27308) @@ -0,0 +1,66 @@ +%Test Name: YiSystemDebug +%Forward Love number solution for a viscoelastic earth, model M3-L70-V01 from +%Spada, G., Barletta, V. R., Klemann, V., Riva, R. E. M., Martinec, Z., +%Gasperini, P., Lund, B., Wolf, D., Vermeersen, L. L. A. and King, M. A. +%(2011), A benchmark study for glacial isostatic adjustment codes. Geophysical +%Journal International, 185: 106--132. doi:10.1111/j.1365-246X.2011.04952.x + +md=model(); +md.cluster=generic('name',oshostname(),'np',2); + +md.materials=materials('litho'); +md.miscellaneous.name='test2073'; +md.groundingline.migration='None'; + +md.verbose=verbose('all'); +md.verbose=verbose('1111111111111111'); +cst=365.25*24*3600*1000; + +md.materials.numlayers=6; +md.materials.radius = [10 1222.5 3.4800e+03 5.7010e+03 5.9510e+03 6.3010e+03 6.3710e+03]'*1e3; +md.materials.density= [1.0750e4 1.0750e+04 4.9780e+03 3.8710e+03 3.4380e+03 3.0370e+03]'; +md.materials.lame_mu= [1e-5 0 2.2834e+00 1.0549e+00 7.0363e-01 5.0605e-01]'*1e11; +md.materials.viscosity=[0 0 2.0000e+00 1.0000e+00 1.0000e+00 1.0000e+25]'*1e21; +md.materials.lame_lambda=md.materials.lame_mu*0+5e17; +md.materials.issolid=[1 0 1 1 1 1]'; +md.materials.rheologymodel=zeros(md.materials.numlayers,1); +md.materials.burgers_mu=md.materials.lame_mu/3; +md.materials.burgers_viscosity=md.materials.viscosity/10; +md.materials.ebm_alpha= ones(md.materials.numlayers,1)*.9; +md.materials.ebm_delta= ones(md.materials.numlayers,1)*0.2; +md.materials.ebm_taul= ones(md.materials.numlayers,1)*54*60; %54min +md.materials.ebm_tauh= ones(md.materials.numlayers,1)*18.6*cst/1e3; %18.6yr + +md.love.allow_layer_deletion=1; +md.love.frequencies=[0]; +md.love.nfreq=length(md.love.frequencies); +md.love.sh_nmin=2; +md.love.sh_nmax=2; +md.love.underflow_tol=1e-20; +md.love.pw_threshold=1e-3; +md.love.Gravitational_Constant=6.6732e-11; +md.love.allow_layer_deletion=1; +md.love.forcing_type=11; +md.love.chandler_wobble=0; +md.love.complex_computation=0; +md.love.istemporal=0; + +md.love.debug=1; +%make sure we only have one step per layer in the propagator +md.love.max_integration_dr=md.materials.radius(end); +md.love.min_integration_steps=1; + +md.love.love_kernels=1; + +md=solve(md,'lv'); + +yi=md.results.LoveSolution.LoveYi; +rhs=md.results.LoveSolution.LoveRhs; + + +%Fields and tolerances to track changes + +field_names ={'LoveYi_loading_elastic','LoveRhs_loading_elastic'}; +field_tolerances={1.0e-13,1.0e-13,1.0e-13}; +field_values={yi,rhs}; + Index: /issm/trunk-jpl/test/NightlyRun/test2084.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2084.m (revision 27307) +++ /issm/trunk-jpl/test/NightlyRun/test2084.m (revision 27308) @@ -44,5 +44,4 @@ md.love.pw_threshold=1e-3; md.love.Gravitational_Constant=6.6732e-11; -md.love.integration_steps_per_layer=100; md.love.allow_layer_deletion=1; md.love.forcing_type=11;