Index: sm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1301.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1301.py (revision 14101) +++ (revision ) @@ -1,87 +1,0 @@ -""" -== == == == == == == == == == == == == == == == == == == -Auto generated python script for ISSM: test1301.m -Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep -== == == == == == == == == == == == == == == == == == == - -Matlab script conversion into python -translateToPy.py Author: Michael Pellegrin -translateToPy.py Date: 09/24/12 -== == == == == == == == == == == == == == == == == == == -""" - -from MatlabFuncs import * -from model import * -from EnumDefinitions import * -from numpy import * -from triangle import * -from setmask import * -from parameterize import * -from setflowequation import * -from solve import * - -# This file can be run to check that the melting in simple conduction is correctly modeled. - -# There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature - -# at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower - -# surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed... - -printingflag=false - - -md=model() -md=triangle(md,'../Exp/Square.exp',100000) -md=setmask(md,'','') -md=parameterize(md,'../Par/SquareThermal.py') -md.extrude(3,2.) -md=setflowequation(md,'Pattyn','all') - - -# Some conditions specific to melting test - -md.initialization.pressure=zeros(md.mesh.numberofvertices,1) -md.initialization.temperature=273.15*ones(md.mesh.numberofvertices,1) -pos=numpy.nonzero(md.mesh.vertexonsurface) -md.thermal.spctemperature(pos)=md.initialization.temperature(pos) -md.materials.rheology_B=paterson(md.initialization.temperature) - - -# analytical results - -# melting heat = geothermal flux - -# Mb*L*rho=G => Mb=G/L*rho - -melting=md.basalforcings.geothermalflux/(md.materials.rho_ice*md.materials.latentheat)*md.constants.yts - - -# modeled results - -md.cluster=generic('name',oshostname(),'np',2) -md=solve(md,ThermalSolutionEnum()) - - -# plot results - -comp_melting=md.results.ThermalSolution.BasalforcingsMeltingRate -relative=abs((comp_melting-melting)./melting)*100 -relative(find(comp_melting==melting))=0 -plotmodel()(md,'data',comp_melting,'title','Modeled melting','data',melting,'title','Analytical melting',\ - 'data',comp_melting-melting,'title','Absolute error','data',relative,'title','Relative error [%]',\ - 'layer#all',1,'caxis#2',[1.02964 1.02966]*10^-4,'FontSize#all',20,'figposition','mathieu') -if printingflag, - set(gcf,'Color','w') - printmodel()('thermalmelting','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off') - system(['mv thermalmelting.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']) -end - - - - -# Fields and tolerances to track changes - -field_names =['BasalMelting'] -field_tolerances=[1e-08] -field_values =[comp_melting] Index: sm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1302.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1302.py (revision 14101) +++ (revision ) @@ -1,86 +1,0 @@ -""" -== == == == == == == == == == == == == == == == == == == -Auto generated python script for ISSM: test1302.m -Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep -== == == == == == == == == == == == == == == == == == == - -Matlab script conversion into python -translateToPy.py Author: Michael Pellegrin -translateToPy.py Date: 09/24/12 -== == == == == == == == == == == == == == == == == == == -""" - -from MatlabFuncs import * -from model import * -from EnumDefinitions import * -from numpy import * -from triangle import * -from setmask import * -from parameterize import * -from setflowequation import * -from solve import * - -# This file can be run to check that the advection-diffusion is correctly modeled. - -# There is u=v=0 and w=cst everywhere the only thermal boundary conditions are an imposed temperature - -# at upper surface and an impose flux at its base. - -printingflag=false - - -md=model() -md=triangle(md,'../Exp/Square.exp',100000) -md=setmask(md,'','') -md=parameterize(md,'../Par/SquareThermal.py') -md.extrude(30,1.) %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%.) -md=setflowequation(md,'Pattyn','all') - - -# Thermal boundary conditions - -pos1=numpy.nonzero(md.mesh.elementonbed) md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10 -pos2=numpy.nonzero(md.mesh.elementonsurface) md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0 -md.initialization.vz=0.1*ones(md.mesh.numberofvertices,1) -md.initialization.vel=sqrt( md.initialization.vx.^2+ md.initialization.vy.^2+ md.initialization.vz.^2) -md.initialization.pressure=zeros(md.mesh.numberofvertices,1) - - -md.thermal.stabilization=2 -# analytical results - -# d2T/dz2-w*rho_ice*c/k*dT/dz=0 T(surface)=0 T(bed)=10 => T=A exp(alpha z)+B - -alpha=0.1/md.constants.yts*md.materials.rho_ice*md.materials.heatcapacity/md.materials.thermalconductivity %alpha=w rho_ice c /k and w=0.1m/an -A=10/(exp(alpha*(-1000))-1) %A=T(bed)/(exp(alpha*bed)-1) with bed=-1000 T(bed)=10 -B=-A -md.initialization.temperature=A*exp(alpha*md.mesh.z)+B - - -# modeled results - -md.cluster=generic('name',oshostname(),'np',2) -md=solve(md,ThermalSolutionEnum()) - - -# plot results - -comp_temp=md.results.ThermalSolution.Temperature -relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100 -relative(find(comp_temp==md.initialization.temperature))=0 -plotmodel()(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,\ - 'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,\ - 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,\ - 'figposition','mathieu','FontSize#all',20) -if printingflag, - set(gcf,'Color','w') - printmodel()('thermaladvection','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off') - system(['mv thermaladvection.png ' ISSM_DIR '/website/doc_pdf/validation/Images/EISMINT ']) -end - - -# Fields and tolerances to track changes - -field_names =['AdvectionTemperature'] -field_tolerances=[1e-13] -field_values =[comp_temp] Index: sm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1303.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1303.py (revision 14101) +++ (revision ) @@ -1,80 +1,0 @@ -""" -== == == == == == == == == == == == == == == == == == == -Auto generated python script for ISSM: test1303.m -Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep -== == == == == == == == == == == == == == == == == == == - -Matlab script conversion into python -translateToPy.py Author: Michael Pellegrin -translateToPy.py Date: 09/24/12 -== == == == == == == == == == == == == == == == == == == -""" - -from MatlabFuncs import * -from model import * -from EnumDefinitions import * -from numpy import * -from triangle import * -from setmask import * -from parameterize import * -from setflowequation import * -from solve import * - -# This file can be run to check that the conduction is correctly modeled. - -# There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature - -# at the lower and upper surface. The result must be a linear temperature from the upper to the lower - -# surface. if it is not the case, something is thermal modeling has been changed... - -printingflag=false - - -md=model() -md=triangle(md,'../Exp/Square.exp',100000) -md=setmask(md,'all','') -md=parameterize(md,'../Par/SquareThermal.py') -md.extrude(11,2.) -md=setflowequation(md,'Pattyn','all') -pos1=numpy.nonzero(md.mesh.elementonbed) md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10 -pos2=numpy.nonzero(md.mesh.elementonsurface) md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0 -md.initialization.pressure=zeros(md.mesh.numberofvertices,1) - - -# analytical results - -# d2T/dz2=0 T(bed)=10 T(surface)=0 => T=0*(z-bed)/thickness+10*(surface-z)/thickness - -# each layer of the 3d mesh must have a constant value - -md.initialization.temperature=10*(md.geometry.surface-md.mesh.z)./md.geometry.thickness - - -# modeled results - -md.cluster=generic('name',oshostname(),'np',2) -md=solve(md,ThermalSolutionEnum()) - - -# plot results - -comp_temp=md.results.ThermalSolution.Temperature -relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100 -relative(find(comp_temp==md.initialization.temperature))=0 -plotmodel()(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,\ - 'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,\ - 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,\ - 'figposition','mathieu','FontSize#all',20) -if printingflag, - set(gcf,'Color','w') - printmodel()('thermalconduction','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off') - system(['mv thermalconduction.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']) -end - - -# Fields and tolerances to track changes - -field_names =['ConductionTemperature'] -field_tolerances=[1e-13] -field_values =[comp_temp] Index: sm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1304.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1304.py (revision 14101) +++ (revision ) @@ -1,82 +1,0 @@ -""" -== == == == == == == == == == == == == == == == == == == -Auto generated python script for ISSM: test1304.m -Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep -== == == == == == == == == == == == == == == == == == == - -Matlab script conversion into python -translateToPy.py Author: Michael Pellegrin -translateToPy.py Date: 09/24/12 -== == == == == == == == == == == == == == == == == == == -""" - -from MatlabFuncs import * -from model import * -from EnumDefinitions import * -from numpy import * -from triangle import * -from setmask import * -from parameterize import * -from setflowequation import * -from solve import * - -# This file can be run to check that the geothermal flux in simple conduction is correctly modeled. - -# There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature - -# at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower - -# surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed... - -printingflag=false - - -md=model() -md=triangle(md,'../Exp/Square.exp',100000) -md=setmask(md,'','') -md=parameterize(md,'../Par/SquareThermal.py') -md.extrude(11,1.) -md=setflowequation(md,'Pattyn','all') - - -pos2=numpy.nonzero(md.mesh.elementonsurface) md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0 -md.initialization.pressure=zeros(md.mesh.numberofvertices,1) -md.basalforcings.geothermalflux(:)=0.1 %100mW/m^2 - - -# analytical results - -# the result is linear with depth and is equal to 0 on the upper surface (See BC) - -# d2T/dz2=0 -k*dT/dz(bed)=G T(surface)=0 => T=-G/k*(z-surface) - -md.initialization.temperature=-0.1/md.materials.thermalconductivity*(md.mesh.z-md.geometry.surface) %G=0.1 W/m2 - - -# modeled results - -md.cluster=generic('name',oshostname(),'np',2) -md=solve(md,ThermalSolutionEnum()) - - -# plot results - -comp_temp=md.results.ThermalSolution.Temperature -relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100 -relative(find(comp_temp==md.initialization.temperature))=0 -plotmodel()(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,\ - 'title','Analytical temperature','view',3,'data',comp_temp-md.initialization.temperature,\ - 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,\ - 'figposition','mathieu','FontSize#all',20) -if printingflag, - set(gcf,'Color','w') - printmodel()('thermalgeothermalflux','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off') - system(['mv thermalgeothermalflux.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']) -end - - -# Fields and tolerances to track changes - -field_names =['GeothermalFluxTemperature'] -field_tolerances=[1e-13] -field_values =[comp_temp] Index: /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt =================================================================== --- /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt (revision 14101) +++ /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt (revision 14102) @@ -11,2 +11,4 @@ test418 needs Dakota test420 needs Dakota +test1401 roundoff error in metric causes different meshes from matlab +test1402 roundoff error in metric causes different meshes from matlab Index: /issm/trunk-jpl/test/NightlyRun/test1301.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1301.m (revision 14101) +++ /issm/trunk-jpl/test/NightlyRun/test1301.m (revision 14102) @@ -10,5 +10,5 @@ md=parameterize(md,'../Par/SquareThermal.par'); md=extrude(md,3,2.); -md=setflowequation(md,'Pattyn','all'); +md=setflowequation(md,'pattyn','all'); %Some conditions specific to melting test @@ -24,5 +24,5 @@ melting=md.basalforcings.geothermalflux/(md.materials.rho_ice*md.materials.latentheat)*md.constants.yts; -%modeled results +%modeled results md.cluster=generic('name',oshostname(),'np',2); md=solve(md,ThermalSolutionEnum()); @@ -30,10 +30,10 @@ %plot results comp_melting=md.results.ThermalSolution.BasalforcingsMeltingRate; -relative=abs((comp_melting-melting)./melting)*100; -relative(find(comp_melting==melting))=0; +relative=abs((comp_melting-melting)./melting)*100.; +relative(find(comp_melting==melting))=0.; plotmodel(md,'data',comp_melting,'title','Modeled melting','data',melting,'title','Analytical melting',... 'data',comp_melting-melting,'title','Absolute error','data',relative,'title','Relative error [%]',... 'layer#all',1,'caxis#2',[1.02964 1.02966]*10^-4,'FontSize#all',20,'figposition','mathieu') -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('thermalmelting','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off'); @@ -41,5 +41,4 @@ end - %Fields and tolerances to track changes field_names ={'BasalMelting'}; Index: /issm/trunk-jpl/test/NightlyRun/test1301.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1301.py (revision 14102) +++ /issm/trunk-jpl/test/NightlyRun/test1301.py (revision 14102) @@ -0,0 +1,60 @@ +import numpy +from model import * +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from paterson import * +from EnumDefinitions import * +from solve import * +from MatlabFuncs import * + +""" +This file can be run to check that the melting in simple conduction is correctly modeled. +There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature +at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower +surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed... +""" + +printingflag=False + +md=model() +md=triangle(md,'../Exp/Square.exp',100000.) +md=setmask(md,'','') +md=parameterize(md,'../Par/SquareThermal.py') +md.extrude(3,2.) +md=setflowequation(md,'pattyn','all') + +#Some conditions specific to melting test +md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1),int) +md.initialization.temperature=273.15*numpy.ones((md.mesh.numberofvertices,1)) +pos=numpy.nonzero(md.mesh.vertexonsurface)[0] +md.thermal.spctemperature[pos]=md.initialization.temperature[pos] +md.materials.rheology_B=paterson(md.initialization.temperature) + +#analytical results +#melting heat = geothermal flux +#Mb*L*rho=G => Mb=G/L*rho +melting=md.basalforcings.geothermalflux/(md.materials.rho_ice*md.materials.latentheat)*md.constants.yts + +#modeled results +md.cluster=generic('name',oshostname(),'np',2) +md=solve(md,ThermalSolutionEnum()) + +#plot results +comp_melting=md.results.ThermalSolution.BasalforcingsMeltingRate +relative=numpy.abs((comp_melting-melting)/melting)*100. +relative[numpy.nonzero(comp_melting==melting)[0]]=0. +#plotmodel(md,'data',comp_melting,'title','Modeled melting','data',melting,'title','Analytical melting',... +# 'data',comp_melting-melting,'title','Absolute error','data',relative,'title','Relative error [%]',... +# 'layer#all',1,'caxis#2',[1.02964 1.02966]*10^-4,'FontSize#all',20,'figposition','mathieu') +if printingflag: + pass +# set(gcf,'Color','w') +# printmodel('thermalmelting','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off'); +# system(['mv thermalmelting.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']); + +#Fields and tolerances to track changes +field_names =['BasalMelting'] +field_tolerances=[1e-08] +field_values =[comp_melting] Index: /issm/trunk-jpl/test/NightlyRun/test1302.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1302.m (revision 14101) +++ /issm/trunk-jpl/test/NightlyRun/test1302.m (revision 14102) @@ -9,9 +9,9 @@ md=parameterize(md,'../Par/SquareThermal.par'); md=extrude(md,30,1.); %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%) -md=setflowequation(md,'Pattyn','all'); +md=setflowequation(md,'pattyn','all'); %Thermal boundary conditions -pos1=find(md.mesh.elementonbed); md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10; -pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0; +pos1=find(md.mesh.elementonbed); md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10.; +pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0.; md.initialization.vz=0.1*ones(md.mesh.numberofvertices,1); md.initialization.vel=sqrt( md.initialization.vx.^2+ md.initialization.vy.^2+ md.initialization.vz.^2); @@ -22,9 +22,9 @@ %d2T/dz2-w*rho_ice*c/k*dT/dz=0 T(surface)=0 T(bed)=10 => T=A exp(alpha z)+B alpha=0.1/md.constants.yts*md.materials.rho_ice*md.materials.heatcapacity/md.materials.thermalconductivity; %alpha=w rho_ice c /k and w=0.1m/an -A=10/(exp(alpha*(-1000))-1); %A=T(bed)/(exp(alpha*bed)-1) with bed=-1000 T(bed)=10 +A=10./(exp(alpha*(-1000.))-1.); %A=T(bed)/(exp(alpha*bed)-1) with bed=-1000 T(bed)=10 B=-A; md.initialization.temperature=A*exp(alpha*md.mesh.z)+B; -%modeled results +%modeled results md.cluster=generic('name',oshostname(),'np',2); md=solve(md,ThermalSolutionEnum()); @@ -32,11 +32,11 @@ %plot results comp_temp=md.results.ThermalSolution.Temperature; -relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100; -relative(find(comp_temp==md.initialization.temperature))=0; +relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100.; +relative(find(comp_temp==md.initialization.temperature))=0.; plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,... 'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,... 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,... 'figposition','mathieu','FontSize#all',20) -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('thermaladvection','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off'); Index: /issm/trunk-jpl/test/NightlyRun/test1302.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1302.py (revision 14102) +++ /issm/trunk-jpl/test/NightlyRun/test1302.py (revision 14102) @@ -0,0 +1,65 @@ +import numpy +import sys +from model import * +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from EnumDefinitions import * +from solve import * +from MatlabFuncs import * + +""" +This file can be run to check that the advection-diffusion is correctly modeled. +There is u=v=0 and w=cst everywhere the only thermal boundary conditions are an imposed temperature +at upper surface and an impose flux at its base. +""" + +printingflag=False + +md=model() +md=triangle(md,'../Exp/Square.exp',100000.) +md=setmask(md,'','') +md=parameterize(md,'../Par/SquareThermal.py') +md.extrude(30,1.) #NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%) +md=setflowequation(md,'pattyn','all') + +#Thermal boundary conditions +pos1=numpy.nonzero(md.mesh.elementonbed)[0] +md.thermal.spctemperature[md.mesh.elements[pos1,0:3]-1]=10. +pos2=numpy.nonzero(md.mesh.elementonsurface)[0] +md.thermal.spctemperature[md.mesh.elements[pos2,3:6]-1]=0. +md.initialization.vz=0.1*numpy.ones((md.mesh.numberofvertices,1)) +md.initialization.vel=numpy.sqrt( md.initialization.vx**2+ md.initialization.vy**2+ md.initialization.vz**2) +md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1),int) + +md.thermal.stabilization=2 +#analytical results +#d2T/dz2-w*rho_ice*c/k*dT/dz=0 T(surface)=0 T(bed)=10 => T=A exp(alpha z)+B +alpha=0.1/md.constants.yts*md.materials.rho_ice*md.materials.heatcapacity/md.materials.thermalconductivity #alpha=w rho_ice c /k and w=0.1m/an +A=10./(numpy.exp(alpha*(-1000.))-1.) #A=T(bed)/(exp(alpha*bed)-1) with bed=-1000 T(bed)=10 +B=-A +md.initialization.temperature=A*numpy.exp(alpha*md.mesh.z)+B + +#modeled results +md.cluster=generic('name',oshostname(),'np',2) +md=solve(md,ThermalSolutionEnum()) + +#plot results +comp_temp=md.results.ThermalSolution.Temperature +relative=numpy.abs((comp_temp-md.initialization.temperature)/md.initialization.temperature)*100. +relative[numpy.nonzero(comp_temp==md.initialization.temperature)[0]]=0. +#plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,... +# 'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,... +# 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,... +# 'figposition','mathieu','FontSize#all',20) +if printingflag: + pass +# set(gcf,'Color','w') +# printmodel('thermaladvection','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off'); +# system(['mv thermaladvection.png ' ISSM_DIR '/website/doc_pdf/validation/Images/EISMINT ']); + +#Fields and tolerances to track changes +field_names =['AdvectionTemperature'] +field_tolerances=[1e-13] +field_values =[comp_temp] Index: /issm/trunk-jpl/test/NightlyRun/test1303.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1303.m (revision 14101) +++ /issm/trunk-jpl/test/NightlyRun/test1303.m (revision 14102) @@ -10,7 +10,8 @@ md=parameterize(md,'../Par/SquareThermal.par'); md=extrude(md,11,2.); -md=setflowequation(md,'Pattyn','all'); -pos1=find(md.mesh.elementonbed); md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10; -pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0; +md=setflowequation(md,'pattyn','all'); + +pos1=find(md.mesh.elementonbed); md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10.; +pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0.; md.initialization.pressure=zeros(md.mesh.numberofvertices,1); @@ -18,7 +19,7 @@ %d2T/dz2=0 T(bed)=10 T(surface)=0 => T=0*(z-bed)/thickness+10*(surface-z)/thickness %each layer of the 3d mesh must have a constant value -md.initialization.temperature=10*(md.geometry.surface-md.mesh.z)./md.geometry.thickness; +md.initialization.temperature=10.*(md.geometry.surface-md.mesh.z)./md.geometry.thickness; -%modeled results +%modeled results md.cluster=generic('name',oshostname(),'np',2); md=solve(md,ThermalSolutionEnum()); @@ -26,11 +27,11 @@ %plot results comp_temp=md.results.ThermalSolution.Temperature; -relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100; -relative(find(comp_temp==md.initialization.temperature))=0; +relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100.; +relative(find(comp_temp==md.initialization.temperature))=0.; plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,... 'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,... 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,... 'figposition','mathieu','FontSize#all',20) -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('thermalconduction','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off'); Index: /issm/trunk-jpl/test/NightlyRun/test1303.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1303.py (revision 14102) +++ /issm/trunk-jpl/test/NightlyRun/test1303.py (revision 14102) @@ -0,0 +1,60 @@ +import numpy +import sys +from model import * +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from EnumDefinitions import * +from solve import * +from MatlabFuncs import * + +""" +This file can be run to check that the conduction is correctly modeled. +There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature +at the lower and upper surface. The result must be a linear temperature from the upper to the lower +surface. if it is not the case, something is thermal modeling has been changed... +""" + +printingflag=False + +md=model() +md=triangle(md,'../Exp/Square.exp',100000.) +md=setmask(md,'all','') +md=parameterize(md,'../Par/SquareThermal.py') +md.extrude(11,2.) +md=setflowequation(md,'pattyn','all') + +pos1=numpy.nonzero(md.mesh.elementonbed)[0] +md.thermal.spctemperature[md.mesh.elements[pos1,0:3]-1]=10. +pos2=numpy.nonzero(md.mesh.elementonsurface)[0] +md.thermal.spctemperature[md.mesh.elements[pos2,3:6]-1]=0. +md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1),int) + +#analytical results +#d2T/dz2=0 T(bed)=10 T(surface)=0 => T=0*(z-bed)/thickness+10*(surface-z)/thickness +#each layer of the 3d mesh must have a constant value +md.initialization.temperature=10.*(md.geometry.surface-md.mesh.z)/md.geometry.thickness + +#modeled results +md.cluster=generic('name',oshostname(),'np',2) +md=solve(md,ThermalSolutionEnum()) + +#plot results +comp_temp=md.results.ThermalSolution.Temperature +relative=numpy.abs((comp_temp-md.initialization.temperature)/md.initialization.temperature)*100. +relative[numpy.nonzero(comp_temp==md.initialization.temperature)[0]]=0. +#plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,... +# 'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,... +# 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,... +# 'figposition','mathieu','FontSize#all',20) +if printingflag: + pass +# set(gcf,'Color','w') +# printmodel('thermalconduction','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off'); +# system(['mv thermalconduction.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']); + +#Fields and tolerances to track changes +field_names =['ConductionTemperature'] +field_tolerances=[1e-13] +field_values =[comp_temp] Index: /issm/trunk-jpl/test/NightlyRun/test1304.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1304.m (revision 14101) +++ /issm/trunk-jpl/test/NightlyRun/test1304.m (revision 14102) @@ -10,7 +10,7 @@ md=parameterize(md,'../Par/SquareThermal.par'); md=extrude(md,11,1.); -md=setflowequation(md,'Pattyn','all'); +md=setflowequation(md,'pattyn','all'); -pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0; +pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0.; md.initialization.pressure=zeros(md.mesh.numberofvertices,1); md.basalforcings.geothermalflux(:)=0.1; %100mW/m^2 @@ -21,5 +21,5 @@ md.initialization.temperature=-0.1/md.materials.thermalconductivity*(md.mesh.z-md.geometry.surface); %G=0.1 W/m2 -%modeled results +%modeled results md.cluster=generic('name',oshostname(),'np',2); md=solve(md,ThermalSolutionEnum()); @@ -27,11 +27,11 @@ %plot results comp_temp=md.results.ThermalSolution.Temperature; -relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100; -relative(find(comp_temp==md.initialization.temperature))=0; +relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100.; +relative(find(comp_temp==md.initialization.temperature))=0.; plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,... 'title','Analytical temperature','view',3,'data',comp_temp-md.initialization.temperature,... 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,... 'figposition','mathieu','FontSize#all',20) -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('thermalgeothermalflux','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off'); Index: /issm/trunk-jpl/test/NightlyRun/test1304.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1304.py (revision 14102) +++ /issm/trunk-jpl/test/NightlyRun/test1304.py (revision 14102) @@ -0,0 +1,59 @@ +import numpy +import sys +from model import * +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from EnumDefinitions import * +from solve import * +from MatlabFuncs import * + +""" +This file can be run to check that the geothermal flux in simple conduction is correctly modeled. +There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature +at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower +surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed... +""" + +printingflag=False + +md=model() +md=triangle(md,'../Exp/Square.exp',100000.) +md=setmask(md,'','') +md=parameterize(md,'../Par/SquareThermal.py') +md.extrude(11,1.) +md=setflowequation(md,'pattyn','all') + +pos2=numpy.nonzero(md.mesh.elementonsurface)[0] +md.thermal.spctemperature[md.mesh.elements[pos2,3:6]-1]=0. +md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1),int) +md.basalforcings.geothermalflux[:]=0.1 #100mW/m^2 + +#analytical results +#the result is linear with depth and is equal to 0 on the upper surface (See BC) +#d2T/dz2=0 -k*dT/dz(bed)=G T(surface)=0 => T=-G/k*(z-surface) +md.initialization.temperature=-0.1/md.materials.thermalconductivity*(md.mesh.z-md.geometry.surface) #G=0.1 W/m2 + +#modeled results +md.cluster=generic('name',oshostname(),'np',2) +md=solve(md,ThermalSolutionEnum()) + +#plot results +comp_temp=md.results.ThermalSolution.Temperature +relative=numpy.abs((comp_temp-md.initialization.temperature)/md.initialization.temperature)*100. +relative[numpy.nonzero(comp_temp==md.initialization.temperature)[0]]=0. +#plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,... +# 'title','Analytical temperature','view',3,'data',comp_temp-md.initialization.temperature,... +# 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,... +# 'figposition','mathieu','FontSize#all',20) +if printingflag: + pass +# set(gcf,'Color','w') +# printmodel('thermalgeothermalflux','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off'); +# system(['mv thermalgeothermalflux.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']); + +#Fields and tolerances to track changes +field_names =['GeothermalFluxTemperature'] +field_tolerances=[1e-13] +field_values =[comp_temp] Index: /issm/trunk-jpl/test/NightlyRun/test1401.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1401.m (revision 14101) +++ /issm/trunk-jpl/test/NightlyRun/test1401.m (revision 14102) @@ -1,8 +1,8 @@ %test the anisotropic mesh adaptation -%function to capture = exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; +%function to capture = exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.; printingflag=false; %create square mesh -L=1; %in m +L=1.; %in m nx=70; %numberof nodes in x direction ny=70; @@ -11,7 +11,7 @@ %mesh adaptation loop YAMS md=squaremesh(md,L,L,nx,ny); -md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; +md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.; plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh1_yams1','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -19,8 +19,8 @@ end -md=YamsCall(md,md.inversion.vel_obs,0.001,0.3,1.3,10^-4); -md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; +md=YamsCall(md,md.inversion.vel_obs,0.001,0.3,1.3,10.^-4); +md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.; plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh1_yams2','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -29,7 +29,7 @@ md=YamsCall(md,md.inversion.vel_obs,0.001,0.3,2.5,0.008); -md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; +md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.; plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh1_yams3','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -41,7 +41,7 @@ %mesh adaptation loop BAMG md=squaremesh(md,L,L,nx,ny); -md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; +md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.; plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh1_bamg1','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -50,8 +50,8 @@ md.private.bamg=NaN; -md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',1.3,'err',10^-4); -md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; +md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',1.3,'err',10.^-4); +md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.; plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh1_bamg2','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -61,7 +61,7 @@ md.private.bamg=NaN; md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',2.5,'err',0.008); -md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; +md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.; plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh1_bamg3','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); Index: /issm/trunk-jpl/test/NightlyRun/test1402.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test1402.m (revision 14101) +++ /issm/trunk-jpl/test/NightlyRun/test1402.m (revision 14102) @@ -3,5 +3,5 @@ %create square mesh -L=1; %in m +L=1.; %in m nx=30; %numberof nodes in x direction ny=30; @@ -10,10 +10,10 @@ %mesh adaptation loop YAMS md=squaremesh(md,L,L,nx,ny); -u=4*md.mesh.x-2; v=4*md.mesh.y-2; -md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... - +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... - +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; +u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.; +md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ... + +tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25)); plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh2_yams1','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -21,11 +21,11 @@ end -md=YamsCall(md,md.inversion.vel_obs,0.005,0.3,2.3,10^-2); -u=4*md.mesh.x-2; v=4*md.mesh.y-2; -md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... - +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... - +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; +md=YamsCall(md,md.inversion.vel_obs,0.005,0.3,2.3,10.^2); +u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.; +md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ... + +tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25)); plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh2_yams2','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -34,10 +34,10 @@ md=YamsCall(md,md.inversion.vel_obs,0.005,0.3,3,0.005); -u=4*md.mesh.x-2; v=4*md.mesh.y-2; -md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... - +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... - +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; +u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.; +md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ... + +tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25)); plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh2_yams3','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -49,10 +49,10 @@ %mesh adaptation loop BAMG md=squaremesh(md,L,L,nx,ny); -u=4*md.mesh.x-2; v=4*md.mesh.y-2; -md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... - +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... - +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; +u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.; +md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ... + +tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25)); plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh2_bamg1','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -61,11 +61,11 @@ md.private.bamg=NaN; -md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',2.3,'err',10^-2); -u=4*md.mesh.x-2; v=4*md.mesh.y-2; -md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... - +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... - +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; +md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',2.3,'err',10.^2); +u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.; +md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ... + +tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25)); plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh2_bamg2','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -75,10 +75,10 @@ md.private.bamg=NaN; md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',3,'err',0.005); -u=4*md.mesh.x-2; v=4*md.mesh.y-2; -md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... - +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... - +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; +u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.; +md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ... + +tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25)); plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh2_bamg3','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); @@ -88,10 +88,10 @@ md.private.bamg=NaN; md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',1.5,'err',0.003,'anisomax',1); -u=4*md.mesh.x-2; v=4*md.mesh.y-2; -md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... - +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... - +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; +u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.; +md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ... + +tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25)); plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); -if printingflag, +if printingflag, set(gcf,'Color','w') printmodel('mesh2_bamgiso','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off'); Index: /issm/trunk-jpl/test/Par/Pig.py =================================================================== --- /issm/trunk-jpl/test/Par/Pig.py (revision 14101) +++ /issm/trunk-jpl/test/Par/Pig.py (revision 14102) @@ -48,5 +48,5 @@ md.prognostic.stabilization=1. md.verbose=verbose(0) -md.settings.waitonlock=30. +md.settings.waitonlock=30 md.timestepping.time_step=1. md.timestepping.final_time=2. Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 14101) +++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py (revision 14102) @@ -51,5 +51,5 @@ md.thermal.stabilization=1. md.verbose=verbose(0) -md.settings.waitonlock=30. +md.settings.waitonlock=30 md.diagnostic.restol=0.05 md.steadystate.reltol=0.05 Index: /issm/trunk-jpl/test/Par/SquareSheetShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareSheetShelf.py (revision 14101) +++ /issm/trunk-jpl/test/Par/SquareSheetShelf.py (revision 14102) @@ -57,5 +57,5 @@ md.thermal.stabilization=1 md.verbose=verbose(0) -md.settings.waitonlock=30. +md.settings.waitonlock=30 md.diagnostic.restol=0.05 md.steadystate.reltol=0.05 Index: /issm/trunk-jpl/test/Par/SquareShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelf.py (revision 14101) +++ /issm/trunk-jpl/test/Par/SquareShelf.py (revision 14102) @@ -70,5 +70,5 @@ md.prognostic.stabilization = 1. md.thermal.stabilization = 1. -md.settings.waitonlock = 30. +md.settings.waitonlock = 30 md.verbose=verbose() md.diagnostic.restol = 0.10 Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelfConstrained.py (revision 14101) +++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.py (revision 14102) @@ -55,5 +55,5 @@ md.thermal.stabilization=1. md.verbose = verbose(0) -md.settings.waitonlock=30. +md.settings.waitonlock=30 md.diagnostic.restol=0.05 md.diagnostic.reltol=0.05 Index: /issm/trunk-jpl/test/Par/SquareThermal.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareThermal.par (revision 14101) +++ /issm/trunk-jpl/test/Par/SquareThermal.par (revision 14102) @@ -4,7 +4,7 @@ disp(' creating thickness'); -h=1000; +h=1000.; md.geometry.thickness=h*ones(md.mesh.numberofvertices,1); -md.geometry.bed=-1000*ones(md.mesh.numberofvertices,1); +md.geometry.bed=-1000.*ones(md.mesh.numberofvertices,1); md.geometry.surface=md.geometry.bed+md.geometry.thickness; @@ -15,21 +15,21 @@ disp(' creating drag'); -md.friction.coefficient=200*ones(md.mesh.numberofvertices,1); %q=1. +md.friction.coefficient=200.*ones(md.mesh.numberofvertices,1); %q=1. %Take care of iceshelves: no basal drag pos=find(md.mask.elementonfloatingice); -md.friction.coefficient(md.mesh.elements(pos,:))=0; +md.friction.coefficient(md.mesh.elements(pos,:))=0.; md.friction.p=ones(md.mesh.numberofelements,1); md.friction.q=ones(md.mesh.numberofelements,1); disp(' creating temperatures'); -md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1); +md.initialization.temperature=(273.-20.)*ones(md.mesh.numberofvertices,1); -disp(' creating flow law paramter'); +disp(' creating flow law parameter'); md.materials.rheology_B=paterson(md.initialization.temperature); -md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); disp(' creating surface mass balance'); md.surfaceforcings.mass_balance=ones(md.mesh.numberofvertices,1)/md.constants.yts; %1m/a -md.basalforcings.melting_rate=0*ones(md.mesh.numberofvertices,1)/md.constants.yts; %1m/a +md.basalforcings.melting_rate=0.*ones(md.mesh.numberofvertices,1)/md.constants.yts; %1m/a %Deal with boundary conditions: @@ -41,3 +41,3 @@ md.thermal.spctemperature(:)=md.initialization.temperature; md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1); -pos=find(md.mask.elementongroundedice);md.basalforcings.geothermalflux(md.mesh.elements(pos,:))=1*10^-3; %1 mW/m^2 +pos=find(md.mask.elementongroundedice);md.basalforcings.geothermalflux(md.mesh.elements(pos,:))=1.*10^-3; %1 mW/m^2 Index: /issm/trunk-jpl/test/Par/SquareThermal.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareThermal.py (revision 14102) +++ /issm/trunk-jpl/test/Par/SquareThermal.py (revision 14102) @@ -0,0 +1,48 @@ +import numpy +from paterson import * +from SetMarineIceSheetBC import * + +#Ok, start defining model parameters here + +md.timestepping.time_step=0 + +print " creating thickness" +h=1000. +md.geometry.thickness=h*numpy.ones((md.mesh.numberofvertices,1)) +md.geometry.bed=-1000.*numpy.ones((md.mesh.numberofvertices,1)) +md.geometry.surface=md.geometry.bed+md.geometry.thickness; + +print " creating velocities" +md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1)) +md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1)) +md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1)) + +print " creating drag" +md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1)) #q=1. +#Take care of iceshelves: no basal drag +pos=numpy.nonzero(md.mask.elementonfloatingice)[0] +md.friction.coefficient[md.mesh.elements[pos,:]-1]=0. +md.friction.p=numpy.ones((md.mesh.numberofelements,1)) +md.friction.q=numpy.ones((md.mesh.numberofelements,1)) + +print " creating temperatures" +md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1)) + +print " creating flow law parameter" +md.materials.rheology_B=paterson(md.initialization.temperature) +md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1)) + +print " creating surface mass balance" +md.surfaceforcings.mass_balance=numpy.ones((md.mesh.numberofvertices,1))/md.constants.yts #1m/a +md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices,1))/md.constants.yts #1m/a + +#Deal with boundary conditions: + +print " boundary conditions for diagnostic model" +md=SetMarineIceSheetBC(md,'../Exp/SquareFront.exp') + +print " boundary conditions for thermal model" +md.thermal.spctemperature[:]=md.initialization.temperature +md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) +pos=numpy.nonzero(md.mask.elementongroundedice)[0] +md.basalforcings.geothermalflux[md.mesh.elements[pos,:]-1]=1.*10**-3 #1 mW/m^2