Context Navigation

source: issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp@ 22961

Last change on this file since 22961 was 22961, checked in by Mathieu Morlighem, 7 years ago
CHG: moving some stuff to analyses
File size: 11.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: