Context Navigation

source: issm/trunk/test/NightlyRun/test1501.py@ 20500

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago
merged trunk-jpl and trunk for revision 20497
File size: 7.6 KB

Rev	Line
[19105]	1	#Test Name: SquareShelfTranSawTooth2d
[14022]	2	import numpy
	3	import sys
	4	from model import *
	5	from triangle import *
	6	from setmask import *
	7	from parameterize import *
	8	from setflowequation import *
	9	from EnumDefinitions import *
	10	from solve import *
	11	from MatlabFuncs import *
	12
	13	printingflag = False
	14
	15	md=triangle(model(),'../Exp/Square.exp',350000.)
	16	md=setmask(md,'all','')
	17	md=parameterize(md,'../Par/SquareShelf.py')
[16137]	18	md=setflowequation(md,'SSA','all')
[14022]	19	md.cluster=generic('name',oshostname(),'np',3)
	20	md.transient.isthermal=False
	21
	22	md.timestepping.time_step=1.
	23	md.settings.output_frequency=1
	24	md.timestepping.final_time=2000.
	25
	26	#Solve for thinning rate -> -1 * surface mass balance
	27	smb= 2.*numpy.ones((md.mesh.numberofvertices,1))
[20500]	28	md.smb.mass_balance= smb
[19105]	29	md.basalforcings.groundedice_melting_rate= smb
[14022]	30
[16137]	31	md=solve(md,MasstransportSolutionEnum())
[14022]	32
	33	for i in xrange(1,11):
[16137]	34	md=solve(md,MasstransportSolutionEnum())
[20500]	35	md.smb.mass_balance= md.smb.mass_balance - ((md.results.MasstransportSolution.Thickness)-md.geometry.thickness)
[14022]	36
	37	#Set up transient
[20500]	38	smb = md.smb.mass_balance
[14022]	39
	40	#tooth= [ [ones(400,1)(smb') - 10.]' [ones(400,1)(smb')]' ];
	41	tooth=numpy.hstack((numpy.tile(smb-10.,(1,400)),numpy.tile(smb,(1,400))))
	42	#smb=[ [ones(399,1)*(smb')]' smb tooth tooth];
	43	smb=numpy.hstack((numpy.tile(smb,(1,399)),smb,tooth,tooth))
	44
[20500]	45	#md.smb.mass_balance= smb;
	46	#md.smb.mass_balance(end+1,:)=[1.:2000.];
	47	md.smb.mass_balance=numpy.vstack((smb,numpy.arange(1,2001)))
[14022]	48
	49	md=solve(md,TransientSolutionEnum())
	50
	51	#Fields and tolerances to track changes
[20500]	52	field_names=['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','SmbMassBalance1', \
	53	'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','SmbMassBalance2', \
	54	'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','SmbMassBalance3', \
	55	'Vx4','Vy4','Vel4','Pressure4','Bed4','Surface4','Thickness4','SmbMassBalance4', \
	56	'Vx5','Vy5','Vel5','Pressure5','Bed5','Surface5','Thickness5','SmbMassBalance5']
[14022]	57	field_tolerances=[1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
	58	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
	59	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
	60	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
	61	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10]
	62	field_values=[\
	63	md.results.TransientSolution[400-1].Vx,\
	64	md.results.TransientSolution[400-1].Vy,\
	65	md.results.TransientSolution[400-1].Vel,\
	66	md.results.TransientSolution[400-1].Pressure,\
[17806]	67	md.results.TransientSolution[400-1].Base,\
[14022]	68	md.results.TransientSolution[400-1].Surface,\
	69	md.results.TransientSolution[400-1].Thickness,\
[20500]	70	md.results.TransientSolution[400-1].SmbMassBalance,\
[14022]	71	md.results.TransientSolution[800-1].Vx,\
	72	md.results.TransientSolution[800-1].Vy,\
	73	md.results.TransientSolution[800-1].Vel,\
	74	md.results.TransientSolution[800-1].Pressure,\
[17806]	75	md.results.TransientSolution[800-1].Base,\
[14022]	76	md.results.TransientSolution[800-1].Surface,\
	77	md.results.TransientSolution[800-1].Thickness,\
[20500]	78	md.results.TransientSolution[800-1].SmbMassBalance,\
[14022]	79	md.results.TransientSolution[1200-1].Vx,\
	80	md.results.TransientSolution[1200-1].Vy,\
	81	md.results.TransientSolution[1200-1].Vel,\
	82	md.results.TransientSolution[1200-1].Pressure,\
[17806]	83	md.results.TransientSolution[1200-1].Base,\
[14022]	84	md.results.TransientSolution[1200-1].Surface,\
	85	md.results.TransientSolution[1200-1].Thickness,\
[20500]	86	md.results.TransientSolution[1200-1].SmbMassBalance,\
[14022]	87	md.results.TransientSolution[1600-1].Vx,\
	88	md.results.TransientSolution[1600-1].Vy,\
	89	md.results.TransientSolution[1600-1].Vel,\
	90	md.results.TransientSolution[1600-1].Pressure,\
[17806]	91	md.results.TransientSolution[1600-1].Base,\
[14022]	92	md.results.TransientSolution[1600-1].Surface,\
	93	md.results.TransientSolution[1600-1].Thickness,\
[20500]	94	md.results.TransientSolution[1600-1].SmbMassBalance,\
[14022]	95	md.results.TransientSolution[2000-1].Vx,\
	96	md.results.TransientSolution[2000-1].Vy,\
	97	md.results.TransientSolution[2000-1].Vel,\
	98	md.results.TransientSolution[2000-1].Pressure,\
[17806]	99	md.results.TransientSolution[2000-1].Base,\
[14022]	100	md.results.TransientSolution[2000-1].Surface,\
	101	md.results.TransientSolution[2000-1].Thickness,\
[20500]	102	md.results.TransientSolution[2000-1].SmbMassBalance,\
[14022]	103	]
	104
	105	if printingflag:
	106	pass
	107
	108	"""
	109	starttime = 360;
	110	endtime = 2000;
	111	res = 40;
	112	ts = [starttime:res:endtime];
	113
	114	index = md.mesh.elements;
	115	x1=md.mesh.x(index(:,1)); x2=md.mesh.x(index(:,2)); x3=md.mesh.x(index(:,3));
	116	y1=md.mesh.y(index(:,1)); y2=md.mesh.y(index(:,2)); y3=md.mesh.y(index(:,3));
	117	areas=(0.5((x2-x1).(y3-y1)-(y2-y1).*(x3-x1)));
	118
	119	thickness = [];
	120	volume = [];
	121	massbal = [];
	122	velocity = [];
	123	for t=starttime:endtime
	124	thickness = [thickness (md.results.TransientSolution(t).Thickness)];
	125	volume = [volume mean(md.results.TransientSolution(t).Thickness.value,2).*areas];
[20500]	126	massbal = [massbal (md.results.TransientSolution(t).SmbMassBalance)];
[14022]	127	velocity = [velocity (md.results.TransientSolution(t).Vel)];
	128	end
	129
	130	figure('Position', [0 0 860 932])
	131
	132	options = plotoptions('data','transient_movie','unit','km');
	133	options = options.list{1};
	134	options = checkplotoptions(md,options);
	135
	136	%loop over the time steps
	137	results=md.results.TransientSolution;
	138	count = 1;
	139	for i=ts
	140
	141	subplot(5,9,[28:31 37:40])
	142	set(gca,'pos',get(gca,'pos')+[-0.08 -0.08 0.07 0.08])
	143	field = 'Thickness';
	144
	145	%process data
	146	[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
	147	[data datatype]=processdata(md,results(i).(field),options);
	148
	149	titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
	150	plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
	151	options=changefieldvalue(options,'title',titlestring);
	152	options=addfielddefault(options,'colorbar',1);
	153	options=changefieldvalue(options,'caxis',[0 max(max(thickness))]);
	154	applyoptions(md,[],options);
	155
	156	subplot(5,9,[33:36 42:45])
	157	set(gca,'pos',get(gca,'pos')+[-0.00 -0.08 0.07 0.08])
	158	field = 'Vel';
	159
	160	%process data
	161	[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
	162	[data datatype]=processdata(md,results(i).(field),options);
	163
	164	titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
	165	plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
	166	options=changefieldvalue(options,'title',titlestring);
	167	options=addfielddefault(options,'colorbar',1);
	168	options=changefieldvalue(options,'caxis',[0 max(max(velocity))]);
	169	applyoptions(md,[],options);
	170
	171	subplot(5,4,1:4)
	172	cla
	173	set(gca,'pos',get(gca,'pos')+[-0.07 0.03 0.12 0.015])
	174	plot(starttime:endtime,mean(massbal),'k','LineWidth', 4)
	175	hold on
	176	ya = ylim;
	177	plot([i i], ya, 'r', 'LineWidth',6)
	178	ylim(ya); xlim([starttime endtime]);
	179	title('Surface Mass Balance','FontSize',14)
	180	ylabel('m/year','FontSize',14)
	181
	182	subplot(5,4,5:8)
	183	cla
	184	set(gca,'pos',get(gca,'pos')+[-0.07 0.015 0.12 0.015])
	185	plot(starttime:endtime,sum(volume)/1000/1000/1000,'LineWidth',4)
	186	hold on
	187	ya = ylim;
	188	plot([i i], ya, 'r', 'LineWidth',6)
	189	ylim(ya); xlim([starttime endtime]);
	190	title('Ice Volume','FontSize',14)
	191	ylabel('km^3','FontSize',14)
	192
	193	subplot(5,4,9:12)
	194	cla
	195	set(gca,'pos',get(gca,'pos')+[-0.07 0 0.12 0.015])
	196	plot(starttime:endtime,mean(velocity)/1000, 'LineWidth', 4)
	197	hold on
	198	ya = ylim;
	199	plot([i i], ya, 'r', 'LineWidth',6)
	200	ylim(ya); xlim([starttime endtime]);
	201	title('Mean Velocity','FontSize', 14)
	202	ylabel('km/year','FontSize', 14)
	203	xlabel('year','FontSize', 14)
	204
	205	set(gcf,'Renderer','zbuffer','color','white'); %fixes a bug on Mac OS X (not needed in future Matlab version)
	206	if i==starttime,
	207	%initialize images and frame
	208	frame=getframe(gcf);
	209	[images,map]=rgb2ind(frame.cdata,256,'nodither');
	210	images(1,1,1,length(ts))=0;
	211	else
	212	frame=getframe(gcf);
	213	images(:,:,1,count) = rgb2ind(frame.cdata,map,'nodither');
	214	end
	215
	216	count = count+1;
	217
	218	end
	219
	220	filename='transawtooth2d.gif';
	221	imwrite(images,map,filename,'DelayTime',1.0,'LoopCount',inf)
	222	"""
	223

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

Download in other formats: