Index: /issm/trunk/TEST =================================================================== --- /issm/trunk/TEST (revision 7810) +++ /issm/trunk/TEST (revision 7810) @@ -0,0 +1,1 @@ +this is a test Index: /issm/trunk/examples/Bumps/Bump1_surface_bed/Square.par =================================================================== --- /issm/trunk/examples/Bumps/Bump1_surface_bed/Square.par (revision 7809) +++ /issm/trunk/examples/Bumps/Bump1_surface_bed/Square.par (revision 7810) @@ -13,13 +13,10 @@ %Solution parameters - %parallelization - md.cluster='none'; - md.np=2; - md.time=1; - md.exclusive=0; + md.cluster=eval(oshostname()); %statics md.lowmem=1; - md.eps_rel=0.01; + md.eps_res=0.05; + md.eps_rel=0.1; md.eps_abs=10; %m/yr md.penalty_offset=4; @@ -36,15 +33,4 @@ md.artificial_diffusivity=1; - %control - md.control_type={'drag'}; %'drag', 'B' - md.nsteps=5; - md.tolx=10^-4; - md.maxiter=20; - md.optscal=10; - md.fit='logarithmic'; %'absolute','relative','logarithmic' - md.meanvel=1000/md.yts; %1000 meters/year - md.epsvel=eps; - - disp(' creating thickness'); hmin=500; @@ -58,5 +44,5 @@ %Add bumps to the ice sheet : pos=find(sqrt((md.x-500000).^2+(md.y-500000).^2)<100000); %for a cicular bump - md.bed(pos)=md.bed(pos)+150*cos(pi/200000*(sqrt((md.y(pos)-500000).^2+(md.x(pos)-500000).^2))); %for a circular bump + md.bed(pos)=md.bed(pos)+30*cos(pi/200000*(sqrt((md.y(pos)-500000).^2+(md.x(pos)-500000).^2))); %for a circular bump md.surface=md.bed+md.thickness; %the surface follows the bed so the bump is transmitted to the surface @@ -68,11 +54,10 @@ disp(' creating drag'); - md.drag_type=2; %0 none 1 plastic 2 viscous - md.drag=200*ones(md.numberofgrids,1); %q=1. + md.drag_coefficient=200*ones(md.numberofgrids,1); %q=1. %Take care of iceshelves: no basal drag pos=find(md.elementoniceshelf); - md.drag(md.elements(pos,:))=0; - md.p=ones(md.numberofelements,1); - md.q=ones(md.numberofelements,1); + md.drag_coefficient(md.elements(pos,:))=0; + md.drag_p=ones(md.numberofelements,1); + md.drag_q=ones(md.numberofelements,1); disp(' creating temperature'); @@ -80,47 +65,14 @@ disp(' creating flow law paramter'); - md.B=paterson(md.observed_temperature); - md.n=3*ones(md.numberofelements,1); + md.rheology_B=paterson(md.observed_temperature); + md.rheology_n=3*ones(md.numberofelements,1); disp(' creating accumulation rates'); - md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a - md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a + md.accumulation_rate=ones(md.numberofgrids,1)/md.yts; %1m/a + md.melting_rate=0*ones(md.numberofgrids,1)/md.yts; %1m/a %Deal with boundary conditions: - disp(' boundary conditions for diagnostic model: '); - %Build gridonicefront, array of boundary grids belonging to the icefront: - gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2); - gridonicefront=double(md.gridonboundary & gridinsideicefront); + disp(' boundary conditions'); + md=SetIceShelfBC(md,'Front.exp'); - md.gridondirichlet_diag=zeros(md.numberofgrids,1); - pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1; - md.dirichletvalues_diag=100*ones(md.numberofgrids,2)*[0,0;0,1]; - - pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2))); - md.segmentonneumann_diag=md.segments(pos,:); - md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure) - - disp(' boundary conditions for prognostic model: '); - md.gridondirichlet_prog=zeros(md.numberofgrids,1); - md.dirichletvalues_prog=zeros(md.numberofgrids,1); - pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2))); - md.segmentonneumann_prog=md.segments(pos,:); - md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1); - md.neumannvalues_prog(:)=NaN; %free radiation - - pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2))); - md.segmentonneumann_prog2=md.segments(pos,:); - md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1); - md.neumannvalues_prog2(:)=NaN; %free radiation - - disp(' boundary conditions for thermal model: '); - md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature - md.dirichletvalues_thermal=md.observed_temperature; - md.geothermalflux=zeros(md.numberofgrids,1); - - -% Some Cielo code, ignore. -if strcmp(md.cluster,'yes') - ServerDisconnect; -end Index: /issm/trunk/examples/Bumps/Bump1_surface_bed/runme.m =================================================================== --- /issm/trunk/examples/Bumps/Bump1_surface_bed/runme.m (revision 7809) +++ /issm/trunk/examples/Bumps/Bump1_surface_bed/runme.m (revision 7810) @@ -1,48 +1,22 @@ -% This file can be run to see the response to a bump on the bed of an ice sheet. The bump is not -% transmitted to the surface of the ice. The results of the three models are saved in order to be -% compared. The geometry is square and the ice sheet is flat. The velocity of the icesheet is -% 100 m/yr in one particular direction on every border. Just run this file in Matlab, i -% with a properly setup Ice code. -% Take care to have the same density and number of layers to compare the results easily. - -% Create macayeal model; - -mdm=model; -mdm=mesh(mdm,'DomainOutline.exp',150000); -mdm=geography(mdm,'',''); -mdm=parameterize(mdm,'Square.par'); -mdm=extrude(mdm,6,3); -mdm=setelementstype(mdm,'macayeal','all'); - -%Compute solution with Ice model -mdm=solve(mdm,'diagnostic','ice'); - -save modelmacayeal mdm - -% Create pattyn model; - -mdp=model; -mdp=mesh(mdp,'DomainOutline.exp',150000); -mdp=geography(mdp,'',''); -mdp=parameterize(mdp,'Square.par'); -mdp=extrude(mdp,6,3); -mdp=setelementstype(mdp,'pattyn','all'); - -%Compute solution with Ice model -mdp=solve(mdp,'diagnostic','ice'); - -save modelpattyn mdp - -% Create stokes model; - -mds=model; -mds=mesh(mds,'DomainOutline.exp',150000); -mds=geography(mds,'',''); -mds=parameterize(mds,'Square.par'); -mds=extrude(mds,6,3); -mds=setelementstype(mds,'pattyn','all','stokes','all'); - -%Compute solution with Ice model -mds=solve(mds,'diagnostic','ice'); - -save modelstokes mds +if 0, + md=mesh(model,'DomainOutline.exp',15000); + md=geography(md,'',''); + md=parameterize(md,'Square.par'); + md=extrude(md,2,3); + md=setelementstype(md,'stokes','all'); + md.cluster.np=14; + md.verbose.convergence=1; + md.spcvelocity(:,5)=1000; + md=solve(md,'analysis_type',DiagnosticSolutionEnum); +else + md.vx=PatchToVec(md.results.DiagnosticSolution.Vx); + md.vy=PatchToVec(md.results.DiagnosticSolution.Vy); + md.accumulation_rate(:)=0; + md.melting_rate(:)=0; + md.thickness(:)=1; + md.bed=md.surface-md.thickness; + md.dt=1; + md=solve(md,'analysis_type',PrognosticSolutionEnum); + beep;pause(0.2);beep; + plotmodel(md,'data',PatchToVec(md.results.PrognosticSolution.Thickness)-md.thickness) +end Index: /issm/trunk/externalpackages/export_fig/ghostscript.m =================================================================== --- /issm/trunk/externalpackages/export_fig/ghostscript.m (revision 7809) +++ /issm/trunk/externalpackages/export_fig/ghostscript.m (revision 7810) @@ -124,4 +124,4 @@ function current_gs_path_str = current_gs_path -current_gs_path_str = 'C:\Program Files\gs\gs8.71\bin\gswin32c.exe'; +current_gs_path_str = 'gs'; return Index: /issm/trunk/src/m/utils/Mesh/FixMesh.m =================================================================== --- /issm/trunk/src/m/utils/Mesh/FixMesh.m (revision 7809) +++ /issm/trunk/src/m/utils/Mesh/FixMesh.m (revision 7810) @@ -15,5 +15,4 @@ y2=y; value2=value; - %First, look for orphan vertices, and take them out. @@ -39,5 +38,4 @@ end - %Check all triangles are well oriented. aires=GetAreas(index2,x2,y2); Index: /issm/trunk/src/mex/Test/Test.cpp =================================================================== --- /issm/trunk/src/mex/Test/Test.cpp (revision 7809) +++ /issm/trunk/src/mex/Test/Test.cpp (revision 7810) @@ -11,4 +11,6 @@ char *name = NULL; bool logical; + double *matrix = NULL; + int numel; printf("Parsing options:\n"); @@ -24,9 +26,13 @@ printf("\n value of \"CELL[0,2]\" is %g\n",test); - options->Get(&logical,"logical"); + options->Get(&logical,"logical",true); printf("\n value of \"logical\" is %s\n",logical?"true":"false"); + + options->Get(&matrix,&numel,"matrix"); + printarray(matrix,numel); delete options; xfree((void**)&name); + xfree((void**)&matrix); }