Index: /issm/trunk-jpl/src/m/classes/clusters/generic.py =================================================================== --- /issm/trunk-jpl/src/m/classes/clusters/generic.py (revision 24568) +++ /issm/trunk-jpl/src/m/classes/clusters/generic.py (revision 24569) @@ -30,6 +30,6 @@ self.executionpath = issmdir() + '/execution' self.valgrind = issmdir() + '/externalpackages/valgrind/install/bin/valgrind' - self.valgrindlib = issmdir() + '/externalpackages/valgrind/install/lib/libmpidebug.so' - self.valgrindsup = issmdir() + '/externalpackages/valgrind/issm.supp' + self.valgrindlib = '/usr/lib/x86_64-linux-gnu/openmpi/lib/libmpi.so.20.10.1' + self.valgrindsup = [issmdir() + '/externalpackages/valgrind/issm.supp'] # add any .supp in list form as needed #use provided options to change fields @@ -103,10 +103,14 @@ #Add --gen -suppressions = all to get suppression lines #fid.write('LD_PRELOAD={} \\\n'.format(self.valgrindlib)) it could be deleted + supstring = '' + for supfile in self.valgrindsup: + supstring += ' --suppressions=' + supfile + if IssmConfig('_HAVE_MPI_')[0]: - fid.write('mpiexec -np {} {} --leak-check=full --suppressions={} {}/{} {} {}/{} {} 2>{}.errlog>{}.outlog '. - format(self.np, self.valgrind, self.valgrindsup, self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) + fid.write('mpiexec -np {} {} --leak-check=full {} {}/{} {} {}/{} {} 2>{}.errlog>{}.outlog '. + format(self.np, self.valgrind, supstring, self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) else: - fid.write('{} --leak-check=full --suppressions={} {}/{} {} {}/{} {} 2>{}.errlog>{}.outlog '. - format(self.valgrind, self.valgrindsup, self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) + fid.write('{} --leak-check=full {} {}/{} {} {}/{} {} 2>{}.errlog>{}.outlog '. + format(self.valgrind, supstring, self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) if not io_gather: #concatenate the output files: @@ -147,6 +151,5 @@ #Add - - gen - suppressions = all to get suppression lines #fid.write('LD_PRELOAD={} \\\n'.format(self.valgrindlib)) - fid.write('mpiexec -np {} {} --leak -check=full --suppressions={} {}/kriging.exe {}/{} {} 2 > {}.errlog > {}.outlog ' .format - (self.np, self.valgrind, self.valgrindsup, self.codepath, self.executionpath, modelname, modelname, modelname, modelname)) + fid.write('mpiexec -np {} {} --leak -check=full --suppressions={} {}/kriging.exe {}/{} {} 2 > {}.errlog > {}.outlog ' .format(self.np, self.valgrind, self.valgrindsup, self.codepath, self.executionpath, modelname, modelname, modelname, modelname)) if not io_gather: #concatenate the output files: fid.write('\ncat {}.outbin. *>{}.outbin'.format(modelname, modelname)) Index: /issm/trunk-jpl/src/m/plot/plot_unit.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 24568) +++ /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 24569) @@ -73,4 +73,5 @@ else: limextent = 0. + if options.exist('clim'): lims = options.getfieldvalue('clim', [np.nanmin(data), np.nanmax(data)]) @@ -78,6 +79,7 @@ lims = options.getfieldvalue('caxis', [np.nanmin(data), np.nanmax(data)]) else: - if np.nanmin(data) == np.nanmax(data): - lims = [np.nanmin(data) * 0.9, np.nanmax(data) * 1.1] + if limextent == 0.: + delta = abs(0.1 * np.nanmin(data)) + lims = [np.nanmin(data) - delta, np.nanmax(data) + delta] elif limextent < 1.0e-12: lims = [np.nanmin(data) - 2 * dataspread, np.nanmax(data) + 2 * dataspread] Index: /issm/trunk-jpl/test/NightlyRun/GetIds.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/GetIds.py (revision 24568) +++ /issm/trunk-jpl/test/NightlyRun/GetIds.py (revision 24569) @@ -44,4 +44,6 @@ if np.array([type(i) == int for i in ids_names]).all(): ids = ids_names + if np.array([type(i) == np.int64 for i in ids_names]).all(): + ids = ids_names elif np.array([type(i) == str for i in ids_names]).all(): ids = np.concatenate([IdFromString(i) for i in ids_names]) Index: /issm/trunk-jpl/test/NightlyRun/test900.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test900.m (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test900.m (revision 24569) @@ -0,0 +1,47 @@ +%Test Name:SquareNoDynUnConfinedHydroDC +md=triangle(model(),'../Exp/Square.exp',100000.); +md=setmask(md,'',''); +%reduced slab (20m long) +md.mesh.x=md.mesh.x/5.0e4; +md.mesh.y=md.mesh.y/5.0e4; +md=parameterize(md,'../Par/SquareNoDyn.par'); +md.cluster=generic('name',oshostname(),'np',1); + +md.transient.ishydrology=1; +md.hydrology=(hydrologydc); +md.hydrology=initialize(md.hydrology,md); + +%Hydro Model Parameters +md.hydrology.isefficientlayer=0; +md.hydrology.sedimentlimit_flag=0; +md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1); +md.hydrology.rel_tol=1.0e-6; +md.hydrology.penalty_lock=0; +md.hydrology.max_iter=200; +md.hydrology.transfer_flag=0; +md.hydrology.unconfined_flag=1; +%Sediment +md.hydrology.sediment_porosity=0.1; +md.hydrology.sediment_thickness=10.0; +md.hydrology.sediment_transmitivity=(1.0e-3*md.hydrology.sediment_thickness)*ones(md.mesh.numberofvertices,1); +%init +md.initialization.sediment_head=-5.0*ones(md.mesh.numberofvertices,1); +%BC +md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1); +pos=find(md.mesh.x==0); +md.hydrology.spcsediment_head(pos)=0.5; + +md.timestepping.time_step=5/md.constants.yts; %5s steppin +md.settings.output_frequency=2; +md.timestepping.final_time=300/md.constants.yts; %500s run + +md=solve(md,'Transient'); + +%fields to track, results can also be found in +%Wang 2009 Fig 6b (jouranl of Hydrology) +field_names={'SedimentWaterHead1',... + 'SedimentWaterHead2'}; +field_tolerances={1e-13,... + 1e-13}; +field_values={md.results.TransientSolution(11).SedimentHead,... + md.results.TransientSolution(31).SedimentHead}; Index: /issm/trunk-jpl/test/NightlyRun/test900.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test900.py (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test900.py (revision 24569) @@ -0,0 +1,53 @@ +#Test Name:SquareNoDynUnConfinedHydroDC +import numpy as np +from model import * +from setmask import * +from triangle import triangle +from parameterize import parameterize +from solve import solve +from socket import gethostname +from generic import generic + +md = triangle(model(), '../Exp/Square.exp', 100000.) +md = setmask(md, '', '') +#reduced square (20m long) +md.mesh.x = md.mesh.x / 5.0e4 +md.mesh.y = md.mesh.y / 5.0e4 +md = parameterize(md, '../Par/SquareNoDyn.py') +md.cluster = generic('name', gethostname(), 'np', 1) + +md.transient.ishydrology = True +md.hydrology = hydrologydc() +md.hydrology = md.hydrology.initialize(md) + + +#Hydro Model Parameters +md.hydrology.isefficientlayer = 0 +md.hydrology.sedimentlimit_flag = 0 +md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices)) +md.hydrology.rel_tol = 1.0e-6 +md.hydrology.penalty_lock = 0 +md.hydrology.max_iter = 200 +md.hydrology.transfer_flag = 0 +md.hydrology.unconfined_flag = 1 +#Sediment +md.hydrology.sediment_porosity = 0.1 +md.hydrology.sediment_thickness = 10.0 +md.hydrology.sediment_transmitivity = (1.0e-3 * md.hydrology.sediment_thickness) * np.ones((md.mesh.numberofvertices)) +#init +md.initialization.sediment_head = -5.0 * np.ones((md.mesh.numberofvertices)) +#BC +md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices)) +md.hydrology.spcsediment_head[np.where(md.mesh.x == 0)] = 0.5 + +md.timestepping.time_step = 5 / md.constants.yts #5s steppin +md.settings.output_frequency = 2 +md.timestepping.final_time = 300 / md.constants.yts #500s run + +md = solve(md, 'Transient') + +#fields to track, results can also be found in +#Wang 2009 Fig 6b (journal of Hydrology) +field_names = ['SedimentWaterHead1', 'SedimentWaterHead2'] +field_tolerances = [1e-13, 1e-13] +field_values = [md.results.TransientSolution[10].SedimentHead, md.results.TransientSolution[30].SedimentHead] Index: /issm/trunk-jpl/test/NightlyRun/test901.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test901.m (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test901.m (revision 24569) @@ -0,0 +1,33 @@ +%Test Name: SquareNoDynHydrologyDCOneLayer +md=triangle(model(),'../Exp/Square.exp',100000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/IceCube.par'); +md.cluster=generic('name',oshostname(),'np',1); + +md.transient.ishydrology=1; +md.hydrology=(hydrologydc); +md.hydrology=initialize(md.hydrology,md); + +md.hydrology.isefficientlayer=0; +md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1); +md.hydrology.sedimentlimit_flag=1; +md.hydrology.sedimentlimit=8000.0; + +md.initialization.sediment_head=0.0*ones(md.mesh.numberofvertices,1); +md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1); +pos=find(md.mesh.y==0); +md.hydrology.spcsediment_head(pos)=0.0; +md.basalforcings.groundedice_melting_rate = 2.0*ones(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate = 0.0*ones(md.mesh.numberofvertices,1); +md.hydrology.sediment_transmitivity=3*ones(md.mesh.numberofvertices,1); +md.timestepping.time_step=0; +md.timestepping.final_time=1.0; +md=solve(md,'Hydrology'); + +%Fields and tolerances to track changes +%you can also compare with an analitic solution, but it is exact +%only if no limits are applied +%analitic=(md.mesh.y.^2-2*md.mesh.y*1.0e6)*(-2.0/(2*md.constants.yts*md.hydrology.sediment_transmitivity)) +field_names ={'SedimentWaterHead','SedimentHeadResidual'}; +field_tolerances={1e-13, 3e-10}; +field_values={md.results.HydrologySolution.SedimentHead,md.results.HydrologySolution.SedimentHeadResidual}; Index: /issm/trunk-jpl/test/NightlyRun/test901.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test901.py (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test901.py (revision 24569) @@ -0,0 +1,42 @@ +#Test Name: SquareNoDynHydrologyDCOneLayer +import numpy as np +from model import * +from socket import gethostname +from triangle import triangle +from setmask import setmask +from parameterize import parameterize +from solve import solve +from generic import generic + +md = triangle(model(), '../Exp/Square.exp', 100000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareNoDyn.py') +md.cluster = generic('name', gethostname(), 'np', 1) + +md.transient.ishydrology = True +md.hydrology = hydrologydc() +md.hydrology = md.hydrology.initialize(md) + +md.hydrology.isefficientlayer = 0 +md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices)) +md.hydrology.sedimentlimit_flag = 1 +md.hydrology.sedimentlimit = 8000.0 + +md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices)) +md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices)) +pos = np.nonzero(md.mesh.y == 0.)[0] +md.hydrology.spcsediment_head[pos] = 0.0 +md.basalforcings.groundedice_melting_rate = 2.0 * np.ones((md.mesh.numberofvertices)) +md.basalforcings.floatingice_melting_rate = 0.0 * np.ones((md.mesh.numberofvertices)) +md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices)) +md.timestepping.time_step = 0 +md.timestepping.final_time = 1.0 +md = solve(md, 'Hydrology') + +#Fields and tolerances to track changes +#you can also compare with an analitic solution, but it is exact +#only if no limits are applied +#analitic=(md.mesh.y**2 - 2 * md.mesh.y * 1.0e6) * (-2.0 / (2 * md.constants.yts * md.hydrology.sediment_transmitivity)) +field_names = ['SedimentWaterHead', 'SedimentHeadResidual'] +field_tolerances = [1e-13, 3e-10] +field_values = [md.results.HydrologySolution.SedimentHead, md.results.HydrologySolution.SedimentHeadResidual] Index: /issm/trunk-jpl/test/NightlyRun/test902.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test902.m (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test902.m (revision 24569) @@ -0,0 +1,75 @@ +%Test Name: SquareSheetHydrologyDCTwoLayers +md=triangle(model(),'../Exp/Square.exp',100000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/IceCube.par'); +md.cluster=generic('name',oshostname(),'np',1); + +md.transient.ishydrology=1; +md.hydrology=(hydrologydc); +md.hydrology=initialize(md.hydrology,md); + +md.hydrology.isefficientlayer=1; +md.hydrology.sedimentlimit_flag=1; +md.hydrology.sedimentlimit=800.0; +md.hydrology.transfer_flag = 0; +md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1); +md.initialization.sediment_head=0.0*ones(md.mesh.numberofvertices,1); +md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1); + +md.basalforcings.groundedice_melting_rate = 2.0*ones(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate = 0.0*ones(md.mesh.numberofvertices,1); +md.hydrology.sediment_transmitivity=3*ones(md.mesh.numberofvertices,1); + +md.initialization.epl_head=0.0*ones(md.mesh.numberofvertices,1); +md.initialization.epl_thickness=1.0*ones(md.mesh.numberofvertices,1); +md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1); +md.hydrology.mask_eplactive_node=0*ones(md.mesh.numberofvertices,1); +md.hydrology.epl_conductivity=30; +md.hydrology.epl_initial_thickness=1; +md.hydrology.epl_colapse_thickness=1.0e-3; +md.hydrology.epl_thick_comp=1; +md.hydrology.epl_max_thickness=1; +md.hydrology.steps_per_step=10; +md.timestepping.time_step=2.0; +md.timestepping.final_time=2.0; + +md=solve(md,'Transient'); + +%re-run with no substeps +mdfine=md; +mdfine.results=struct(); +mdfine.hydrology.steps_per_step=1; +mdfine.timestepping.time_step=0.2; +mdfine=solve(mdfine,'Transient'); + +%Fields and tolerances to track changes +field_names ={'SedimentWaterHead1','EplWaterHead1','SedimentHeadResidual1',... + 'SedimentWaterHead4','EplWaterHead4','SedimentHeadResidual4',... + 'SedimentWaterHead5','EplWaterHead5','SedimentHeadResidual5',... + 'SedimentWaterHead9','EplWaterHead9','SedimentHeadResidual9',... + 'EplWaterHead10', 'EplWaterHeadSubstep10', 'SedimentWaterHead10',... + 'SedimentWaterHeadSubstep10'}; +field_tolerances={... + 1e-13, 1e-13, 1e-13,... + 1e-13, 1e-13, 1e-13,... + 1e-13, 5e-12, 1e-11,... + 1e-13, 5e-12, 1e-11,... + 1e-13, 1e-13, 1e-13,... + 1e-13}; +field_values={mdfine.results.TransientSolution(1).SedimentHead, ... + mdfine.results.TransientSolution(1).EplHead,... + mdfine.results.TransientSolution(1).SedimentHeadResidual,... + mdfine.results.TransientSolution(4).SedimentHead,... + mdfine.results.TransientSolution(4).EplHead,... + mdfine.results.TransientSolution(4).SedimentHeadResidual, ... + mdfine.results.TransientSolution(5).SedimentHead,... + mdfine.results.TransientSolution(5).EplHead,... + mdfine.results.TransientSolution(5).SedimentHeadResidual, ... + mdfine.results.TransientSolution(9).SedimentHead,... + mdfine.results.TransientSolution(9).EplHead,... + mdfine.results.TransientSolution(9).SedimentHeadResidual,... + md.results.TransientSolution(1).EplHead,... + md.results.TransientSolution(1).EplHeadSubstep,... + md.results.TransientSolution(1).SedimentHead,... + md.results.TransientSolution(1).SedimentHeadSubstep + }; Index: /issm/trunk-jpl/test/NightlyRun/test902.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test902.py (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test902.py (revision 24569) @@ -0,0 +1,81 @@ +#Test Name: SquareNoDynHydrolgyDCTwoLayers +import numpy as np +from model import * +from socket import gethostname +from triangle import triangle +from setmask import setmask +from parameterize import parameterize +from solve import solve +from generic import generic + +md = triangle(model(), '../Exp/Square.exp', 100000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareNoDyn.py') +md.cluster = generic('name', gethostname(), 'np', 1) + +md.transient.ishydrology = True +md.hydrology = hydrologydc() +md.hydrology = md.hydrology.initialize(md) + +md.hydrology.isefficientlayer = 1 +md.hydrology.sedimentlimit_flag = 1 +md.hydrology.sedimentlimit = 800.0 +md.hydrology.transfer_flag = 0 +md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices)) +md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices)) +md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices)) + +md.basalforcings.groundedice_melting_rate = 2.0 * np.ones((md.mesh.numberofvertices)) +md.basalforcings.floatingice_melting_rate = 0.0 * np.ones((md.mesh.numberofvertices)) +md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices)) + +md.initialization.epl_head = np.zeros((md.mesh.numberofvertices)) +md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices)) +md.hydrology.spcepl_head = np.nan * np.ones((md.mesh.numberofvertices)) +md.hydrology.mask_eplactive_node = np.zeros((md.mesh.numberofvertices)) +md.hydrology.epl_conductivity = 30 +md.hydrology.epl_initial_thickness = 1 +md.hydrology.epl_colapse_thickness = 1.0e-3 +md.hydrology.epl_thick_comp = 1 +md.hydrology.epl_max_thickness = 1 +md.hydrology.steps_per_step = 10 +md.timestepping.time_step = 2.0 +md.timestepping.final_time = 2.0 + +md = solve(md, 'Transient') + +#re-run with no substeps +mdfine = copy.deepcopy(md) +mdfine.results = [] +mdfine.hydrology.steps_per_step = 1 +mdfine.timestepping.time_step = 0.2 +mdfine = solve(mdfine, 'Transient') + +field_names = ['SedimentWaterHead1', 'EplWaterHead1', 'SedimentHeadResidual1', + 'SedimentWaterHead4', 'EplWaterHead4', 'SedimentHeadResidual4', + 'SedimentWaterHead5', 'EplWaterHead5', 'SedimentHeadResidual5', + 'SedimentWaterHead9', 'EplWaterHead9', 'SedimentHeadResidual9', + 'EplWaterHead10', 'EplWaterHeadSubstep10', 'SedimentWaterHead10', + 'SedimentWaterHeadSubstep10'] +field_tolerances = [1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, + 1e-13, 5e-12, 1e-11, + 1e-13, 5e-12, 1e-11, + 1e-13, 1e-13, 1e-13, + 1e-13] +field_values = [mdfine.results.TransientSolution[0].SedimentHead, + mdfine.results.TransientSolution[0].EplHead, + mdfine.results.TransientSolution[0].SedimentHeadResidual, + mdfine.results.TransientSolution[3].SedimentHead, + mdfine.results.TransientSolution[3].EplHead, + mdfine.results.TransientSolution[3].SedimentHeadResidual, + mdfine.results.TransientSolution[4].SedimentHead, + mdfine.results.TransientSolution[4].EplHead, + mdfine.results.TransientSolution[4].SedimentHeadResidual, + mdfine.results.TransientSolution[8].SedimentHead, + mdfine.results.TransientSolution[8].EplHead, + mdfine.results.TransientSolution[8].SedimentHeadResidual, + md.results.TransientSolution[-1].EplHead, + md.results.TransientSolution[-1].EplHeadSubstep, + md.results.TransientSolution[-1].SedimentHead, + md.results.TransientSolution[-1].SedimentHeadSubstep] Index: /issm/trunk-jpl/test/NightlyRun/test903.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test903.m (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test903.m (revision 24569) @@ -0,0 +1,35 @@ +%Test Name: SquareNoDynExtrudedHydrologyDCOneLayer +md=triangle(model(),'../Exp/Square.exp',100000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/IceCube.par'); +md.cluster=generic('name',oshostname(),'np',1); + +md.transient.ishydrology=1; +md.hydrology=(hydrologydc); +md.hydrology=initialize(md.hydrology,md); + +md.hydrology.isefficientlayer=0; +md.hydrology.sedimentlimit_flag=1; +md.hydrology.sedimentlimit=8000.0; +md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1); +md.initialization.sediment_head=0.0*ones(md.mesh.numberofvertices,1); +md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1); +pos=find(md.mesh.y==0); +md.hydrology.spcsediment_head(pos)=0.0; + +md.basalforcings.groundedice_melting_rate = 2.0*ones(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate = 0.0*ones(md.mesh.numberofvertices,1); +md.hydrology.sediment_transmitivity= 3.0*ones(md.mesh.numberofvertices,1); + +md.timestepping.time_step=0; +md.timestepping.final_time=1.0; +md=extrude(md,3,1.1); +md=solve(md,'Hydrology'); + +%Fields and tolerances to track changes +%you can also compare with an analitic solution, but it is exact +%only if no limits are applied +%analitic=(md.mesh.y.^2-2*md.mesh.y*1.0e6)*(-2.0/(2*md.constants.yts*md.hydrology.sediment_transmitivity)) +field_names ={'SedimentWaterHead','SedimentHeadResidual'}; +field_tolerances={1e-13, 3e-10}; +field_values={md.results.HydrologySolution.SedimentHead,md.results.HydrologySolution.SedimentHeadResidual}; Index: /issm/trunk-jpl/test/NightlyRun/test903.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test903.py (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test903.py (revision 24569) @@ -0,0 +1,43 @@ +#Test Name: SquareNoDynExtrudedHydrologyDCOneLayer +import numpy as np +from model import * +from socket import gethostname +from triangle import triangle +from setmask import setmask +from parameterize import parameterize +from solve import solve +from generic import generic + +md = triangle(model(), '../Exp/Square.exp', 100000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareNoDyn.py') +md.cluster = generic('name', gethostname(), 'np', 1) + +md.transient.ishydrology = True +md.hydrology = hydrologydc() +md.hydrology = md.hydrology.initialize(md) + +md.hydrology.isefficientlayer = 0 +md.hydrology.sedimentlimit_flag = 1 +md.hydrology.sedimentlimit = 8000.0 +md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices)) +md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices)) +md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices)) +md.hydrology.spcsediment_head[np.where(md.mesh.y == 0)] = 0.0 + +md.basalforcings.groundedice_melting_rate = 2.0 * np.ones((md.mesh.numberofvertices)) +md.basalforcings.floatingice_melting_rate = 0.0 * np.ones((md.mesh.numberofvertices)) +md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices)) + +md.timestepping.time_step = 0 +md.timestepping.final_time = 1.0 +md.extrude(3, 1.) +md = solve(md, 'Hydrology') + +#Fields and tolerances to track changes +#you can also compare with an analitic solution, but it is exact +#only if no limits are applied +#analitic=(md.mesh.y.^2 - 2 * md.mesh.y * 1.0e6) * (-2.0 / (2 * md.constants.yts * md.hydrology.sediment_transmitivity)) +field_names = ['SedimentWaterHead', 'SedimentHeadResidual'] +field_tolerances = [1e-13, 3e-10] +field_values = [md.results.HydrologySolution.SedimentHead, md.results.HydrologySolution.SedimentHeadResidual] Index: /issm/trunk-jpl/test/NightlyRun/test904.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test904.m (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test904.m (revision 24569) @@ -0,0 +1,57 @@ +%Test Name: SquareNoDynExtrudedHydrologyDCTwoLayers +md=triangle(model(),'../Exp/Square.exp',100000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/IceCube.par'); +md.cluster=generic('name',oshostname(),'np',1); + +md.transient.ishydrology=1; +md.hydrology=(hydrologydc); +md.hydrology=initialize(md.hydrology,md); + +md.hydrology.isefficientlayer=1; +md.hydrology.sedimentlimit_flag=1; +md.hydrology.transfer_flag = 0; +md.hydrology.sedimentlimit=800.0; +md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1); +md.initialization.sediment_head=0.0*ones(md.mesh.numberofvertices,1); +md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1); +md.basalforcings.groundedice_melting_rate = 2.0*ones(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate = 0.0*ones(md.mesh.numberofvertices,1); +md.hydrology.sediment_transmitivity=3*ones(md.mesh.numberofvertices,1); + +md.initialization.epl_head=0.0*ones(md.mesh.numberofvertices,1); +md.initialization.epl_thickness=1.0*ones(md.mesh.numberofvertices,1); +md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1); +md.hydrology.mask_eplactive_node=0*ones(md.mesh.numberofvertices,1); +md.hydrology.epl_conductivity=30; +md.hydrology.epl_initial_thickness=1; +md.hydrology.epl_colapse_thickness=1.0e-3; +md.hydrology.epl_thick_comp=1; +md.hydrology.epl_max_thickness=1; +md.timestepping.time_step=0.2; +md.timestepping.final_time=2.0; + +md=extrude(md,3,1.); +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={'SedimentWaterHead1','EplWaterHead1','SedimentHeadResidual1',... + 'SedimentWaterHead4','EplWaterHead4','SedimentHeadResidual4',... + 'SedimentWaterHead5','EplWaterHead5','SedimentHeadResidual5',... + 'SedimentWaterHead9','EplWaterHead9','SedimentHeadResidual9'}; +field_tolerances={1e-13, 1e-13, 1e-13,... + 1e-13, 1e-13, 1e-13,... + 1e-13, 5e-12, 2e-11,... + 1e-13, 5e-12, 2e-11}; +field_values={md.results.TransientSolution(1).SedimentHead, ... + md.results.TransientSolution(1).EplHead,... + md.results.TransientSolution(1).SedimentHeadResidual,... + md.results.TransientSolution(4).SedimentHead,... + md.results.TransientSolution(4).EplHead,... + md.results.TransientSolution(4).SedimentHeadResidual, ... + md.results.TransientSolution(5).SedimentHead,... + md.results.TransientSolution(5).EplHead,... + md.results.TransientSolution(5).SedimentHeadResidual, ... + md.results.TransientSolution(9).SedimentHead,... + md.results.TransientSolution(9).EplHead,... + md.results.TransientSolution(9).SedimentHeadResidual}; Index: /issm/trunk-jpl/test/NightlyRun/test904.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test904.py (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test904.py (revision 24569) @@ -0,0 +1,66 @@ +#Test Name: SquareNoDynExtrudedHydrologyDCTwoLayers +import numpy as np +from model import * +from socket import gethostname +from triangle import triangle +from setmask import setmask +from parameterize import parameterize +from solve import solve +from generic import generic + +md = triangle(model(), '../Exp/Square.exp', 100000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareNoDyn.py') +md.cluster = generic('name', gethostname(), 'np', 1) + +md.transient.ishydrology = True +md.hydrology = hydrologydc() +md.hydrology = md.hydrology.initialize(md) + +md.hydrology.isefficientlayer = 1 +md.hydrology.sedimentlimit_flag = 1 +md.hydrology.transfer_flag = 0 +md.hydrology.sedimentlimit = 800.0 +md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices)) +md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices)) +md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices)) +md.basalforcings.groundedice_melting_rate = 2.0 * np.ones((md.mesh.numberofvertices)) +md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) +md.hydrology.sediment_transmitivity = 3 * np.ones((md.mesh.numberofvertices)) + +md.initialization.epl_head = np.zeros((md.mesh.numberofvertices)) +md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices)) +md.hydrology.spcepl_head = np.nan * np.ones((md.mesh.numberofvertices)) +md.hydrology.mask_eplactive_node = np.zeros((md.mesh.numberofvertices)) +md.hydrology.epl_conductivity = 30 +md.hydrology.epl_initial_thickness = 1 +md.hydrology.epl_colapse_thickness = 1.0e-3 +md.hydrology.epl_thick_comp = 1 +md.hydrology.epl_max_thickness = 1 +md.timestepping.time_step = 0.2 +md.timestepping.final_time = 2.0 + +md.extrude(3, 1.) +md = solve(md, 'Transient') + +#Fields and tolerances to track changes +field_names = ['SedimentWaterHead1', 'EplWaterHead1', 'SedimentHeadResidual1', + 'SedimentWaterHead4', 'EplWaterHead4', 'SedimentHeadResidual4', + 'SedimentWaterHead5', 'EplWaterHead5', 'SedimentHeadResidual5', + 'SedimentWaterHead9', 'EplWaterHead9', 'SedimentHeadResidual9'] +field_tolerances = [1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, + 1e-13, 5e-12, 2e-11, + 1e-13, 5e-12, 2e-11] +field_values = [md.results.TransientSolution[0].SedimentHead, + md.results.TransientSolution[0].EplHead, + md.results.TransientSolution[0].SedimentHeadResidual, + md.results.TransientSolution[3].SedimentHead, + md.results.TransientSolution[3].EplHead, + md.results.TransientSolution[3].SedimentHeadResidual, + md.results.TransientSolution[4].SedimentHead, + md.results.TransientSolution[4].EplHead, + md.results.TransientSolution[4].SedimentHeadResidual, + md.results.TransientSolution[8].SedimentHead, + md.results.TransientSolution[8].EplHead, + md.results.TransientSolution[8].SedimentHeadResidual] Index: /issm/trunk-jpl/test/NightlyRun/test905.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test905.m (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test905.m (revision 24569) @@ -0,0 +1,84 @@ +%Test Name: SquareNoDynHydrologyDCSmbCoupled +md=triangle(model(),'../Exp/Square.exp',100000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/SquareNoDyn.par'); +md.cluster = generic('name',oshostname(),'np',1); + +md.transient.ishydrology=1; +md.transient.issmb=1; +md.hydrology=(hydrologydc); +md.hydrology=initialize(md.hydrology,md); +md.smb=(SMBgradientscomponents); + +md.hydrology.isefficientlayer=1; + +md.hydrology.sedimentlimit_flag=1; +md.hydrology.sedimentlimit=400.0; +md.hydrology.sediment_transmitivity=3.0*ones(md.mesh.numberofvertices,1); +md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1); + +md.hydrology.mask_eplactive_node=zeros(md.mesh.numberofvertices,1); +md.hydrology.epl_conductivity=3.; +md.hydrology.epl_initial_thickness=20; +md.hydrology.epl_colapse_thickness=1.0e-3; +md.hydrology.epl_thick_comp=0; +md.hydrology.epl_max_thickness=1; + +md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1); +md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1); + +md.initialization.sediment_head=zeros(md.mesh.numberofvertices,1); +md.initialization.epl_head=zeros(md.mesh.numberofvertices,1); +md.initialization.epl_thickness=np.ones(md.mesh.numberofvertices,1); + +md.hydrology.steps_per_step=5; +md.smb.steps_per_step=7; +md.timestepping.time_step=2; +md.timestepping.final_time=20.0; + +smb_step=md.timestepping.time_step/md.smb.steps_per_step; +duration=[md.timestepping.start_time:smb_step:md.timestepping.final_time]; + +ddf=10.0e-3; +md.smb.accuref=[[0.5 0.5];[0. 2.0]]; +md.smb.accualti=0.0; +md.smb.accugrad=0.0; + +md.smb.runoffref=0.9*duration*ddf; +md.smb.runoffref=[md.smb.runoffref;duration]; +md.smb.runoffalti=0.0; +md.smb.runoffgrad=-6.5e-3*ddf; %lapse rate *ddf*day per year + +md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); + +md=solve(md,'Transient'); + +field_names={'SedimentWaterHead1','EplWaterHead1','SedimentHeadResidual1',... + 'SedimentWaterHead4','EplWaterHead4','SedimentHeadResidual4',... + 'SedimentWaterHead5','EplWaterHead5','SedimentHeadResidual5',... + 'SedimentWaterHead9','EplWaterHead9','SedimentHeadResidual9',... + 'EplWaterHead10','EplWaterHeadSubstep10','SedimentWaterHead10',... + 'SedimentWaterHeadSubstep10'} +field_tolerances={1e-13,1e-13,1e-13,... + 1e-13,1e-13,1e-13,... + 1e-13,5e-12,1e-11,... + 1e-13,5e-12,1e-11,... + 1e-13,1e-13,1e-13,... + 1e-13}; +field_values={md.results.TransientSolution(1).SedimentHead,... + md.results.TransientSolution(1).EplHead,... + md.results.TransientSolution(1).SedimentHeadResidual,... + md.results.TransientSolution(4).SedimentHead,... + md.results.TransientSolution(4).EplHead,... + md.results.TransientSolution(4).SedimentHeadResidual,... + md.results.TransientSolution(5).SedimentHead,... + md.results.TransientSolution(5).EplHead,... + md.results.TransientSolution(5).SedimentHeadResidual,... + md.results.TransientSolution(9).SedimentHead,... + md.results.TransientSolution(9).EplHead,... + md.results.TransientSolution(9).SedimentHeadResidual,... + md.results.TransientSolution(end).EplHead,... + md.results.TransientSolution(end).EplHeadSubstep,... + md.results.TransientSolution(end).SedimentHead,... + md.results.TransientSolution(end).SedimentHeadSubstep}; Index: /issm/trunk-jpl/test/NightlyRun/test905.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test905.py (revision 24569) +++ /issm/trunk-jpl/test/NightlyRun/test905.py (revision 24569) @@ -0,0 +1,94 @@ +#Test Name: SquareNoDynHydrologyDCSmbCoupled +import numpy as np +from model import * +from socket import gethostname +from triangle import triangle +from setmask import setmask +from parameterize import parameterize +from solve import solve +from generic import generic +from SMBgradientscomponents import SMBgradientscomponents + +md = triangle(model(), '../Exp/Square.exp', 100000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareNoDyn.py') +md.cluster = generic('name', gethostname(), 'np', 1) + +md.transient.ishydrology = True +md.transient.issmb = True +md.hydrology = hydrologydc() +md.hydrology = md.hydrology.initialize(md) +md.smb = SMBgradientscomponents() + +md.hydrology.isefficientlayer = 1 + +md.hydrology.sedimentlimit_flag = 1 +md.hydrology.sedimentlimit = 400.0 +md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices)) +md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices)) + +md.hydrology.mask_eplactive_node = np.zeros((md.mesh.numberofvertices)) +md.hydrology.epl_conductivity = 3. +md.hydrology.epl_initial_thickness = 20 +md.hydrology.epl_colapse_thickness = 1.0e-3 +md.hydrology.epl_thick_comp = 0 +md.hydrology.epl_max_thickness = 1 + +md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices)) +md.hydrology.spcepl_head = np.nan * np.ones((md.mesh.numberofvertices)) + +md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices)) +md.initialization.epl_head = np.zeros((md.mesh.numberofvertices)) +md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices)) + +md.hydrology.steps_per_step = 5 +md.smb.steps_per_step = 7 #md.hydrology.steps_per_step +md.timestepping.time_step = 2 +md.timestepping.final_time = 20.0 + +smb_step = md.timestepping.time_step / md.smb.steps_per_step +duration = np.arange(md.timestepping.start_time, md.timestepping.final_time + smb_step, smb_step) + +ddf = 10.0e-3 +md.smb.accuref = np.array([[0.5, 0.5], [0., 2.0]]) +md.smb.accualti = 0.0 +md.smb.accugrad = 0.0 + +md.smb.runoffref = 0.9 * duration * ddf +md.smb.runoffref = np.vstack((md.smb.runoffref, duration)) +md.smb.runoffalti = 0.0 +md.smb.runoffgrad = -6.5e-3 * ddf #lapse rate *ddf*day per year + +md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices)) +md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) + +md = solve(md, 'Transient') + +field_names = ['SedimentWaterHead1', 'EplWaterHead1', 'SedimentHeadResidual1', + 'SedimentWaterHead4', 'EplWaterHead4', 'SedimentHeadResidual4', + 'SedimentWaterHead5', 'EplWaterHead5', 'SedimentHeadResidual5', + 'SedimentWaterHead9', 'EplWaterHead9', 'SedimentHeadResidual9', + 'EplWaterHead10', 'EplWaterHeadSubstep10', 'SedimentWaterHead10', + 'SedimentWaterHeadSubstep10'] +field_tolerances = [1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, + 1e-13, 5e-12, 1e-11, + 1e-13, 5e-12, 1e-11, + 1e-13, 1e-13, 1e-13, + 1e-13] +field_values = [md.results.TransientSolution[0].SedimentHead, + md.results.TransientSolution[0].EplHead, + md.results.TransientSolution[0].SedimentHeadResidual, + md.results.TransientSolution[3].SedimentHead, + md.results.TransientSolution[3].EplHead, + md.results.TransientSolution[3].SedimentHeadResidual, + md.results.TransientSolution[4].SedimentHead, + md.results.TransientSolution[4].EplHead, + md.results.TransientSolution[4].SedimentHeadResidual, + md.results.TransientSolution[8].SedimentHead, + md.results.TransientSolution[8].EplHead, + md.results.TransientSolution[8].SedimentHeadResidual, + md.results.TransientSolution[-1].EplHead, + md.results.TransientSolution[-1].EplHeadSubstep, + md.results.TransientSolution[-1].SedimentHead, + md.results.TransientSolution[-1].SedimentHeadSubstep] Index: /issm/trunk-jpl/test/Par/SquareNoDyn.par =================================================================== --- /issm/trunk-jpl/test/Par/SquareNoDyn.par (revision 24569) +++ /issm/trunk-jpl/test/Par/SquareNoDyn.par (revision 24569) @@ -0,0 +1,47 @@ +%Start defining model parameters here + +%Geometry +md.geometry.thickness=1000.0*ones(md.mesh.numberofvertices,1); +md.geometry.base=zeros(md.mesh.numberofvertices,1); +md.geometry.surface=md.geometry.base+md.geometry.thickness; + +%Materials +md.initialization.temperature=(273.-20.)*ones(md.mesh.numberofvertices,1); +md.materials.rheology_B=paterson(md.initialization.temperature); +md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); + +%Friction +md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); +md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.; +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=ones(md.mesh.numberofelements,1); + +%Some necessary fields to fool checkonsistency +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.vz=zeros(md.mesh.numberofvertices,1); +md.initialization.pressure=zeros(md.mesh.numberofvertices,1); + +md.stressbalance.spcvx=zeros(md.mesh.numberofvertices,1); +md.stressbalance.spcvy=zeros(md.mesh.numberofvertices,1); +md.stressbalance.spcvz=zeros(md.mesh.numberofvertices,1); + +md.stressbalance.referential=NaN(md.mesh.numberofvertices,6); +md.stressbalance.loadingforce=zeros(md.mesh.numberofvertices,3); + +md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); + +md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); + + +%Numerical parameters +md.verbose=verbose(0); +md.settings.waitonlock=30; +md.groundingline.migration='None'; + +md.transient=deactivateall(md.transient); + +%Change name so that no test have the same name +A=dbstack; +if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end Index: /issm/trunk-jpl/test/Par/SquareNoDyn.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareNoDyn.py (revision 24569) +++ /issm/trunk-jpl/test/Par/SquareNoDyn.py (revision 24569) @@ -0,0 +1,54 @@ +import os.path +import numpy as np +import inspect +from verbose import verbose +from transient import transient +from paterson import paterson +from arch import * + +#Start defining model parameters here + +#Geometry +md.geometry.thickness = 1000.0 * np.ones((md.mesh.numberofvertices)) +md.geometry.base = np.zeros((md.mesh.numberofvertices)) +md.geometry.surface = md.geometry.base + md.geometry.thickness + +# #Materials +md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices)) +md.materials.rheology_B = paterson(md.initialization.temperature) +md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements)) + +#Friction +md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices)) +md.friction.coefficient[np.where(md.mask.groundedice_levelset < 0.)[0]] = 0. +md.friction.p = np.ones((md.mesh.numberofelements)) +md.friction.q = np.ones((md.mesh.numberofelements)) + +#Some necessary fields to fool checkonsistency +md.initialization.vx = np.zeros((md.mesh.numberofvertices)) +md.initialization.vy = np.zeros((md.mesh.numberofvertices)) +md.initialization.vz = np.zeros((md.mesh.numberofvertices)) +md.initialization.pressure = np.zeros((md.mesh.numberofvertices)) + +md.stressbalance.spcvx = np.zeros((md.mesh.numberofvertices)) +md.stressbalance.spcvy = np.zeros((md.mesh.numberofvertices)) +md.stressbalance.spcvz = np.zeros((md.mesh.numberofvertices)) + +md.stressbalance.referential = float('nan') * np.ones((md.mesh.numberofvertices, 6)) +md.stressbalance.loadingforce = np.zeros((md.mesh.numberofvertices, 3)) + +md.smb.mass_balance = np.zeros((md.mesh.numberofvertices)) + +md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices)) +md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) + +#Numerical parameters +md.verbose = verbose(0) +md.settings.waitonlock = 30 +md.groundingline.migration = 'None' + +md.transient = transient.setallnullparameters(md.transient) + +#Change name so that no test have the same name +if len(inspect.stack()) > 2: + md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]