Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26664) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26665) @@ -1,3 +1,3 @@ -%Test Name: PigTranFrontalforcingsrignotautoregression +%Test Name: PigTranStochasticforcing md=triangle(model(),'../Exp/Pig.exp',10000.); md=setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp'); @@ -6,26 +6,46 @@ md.timestepping.start_time = 0; md.timestepping.time_step = 1; -md.timestepping.final_time = 25; +md.timestepping.final_time = 10; -%separate domain in 2 basins -idbasin = zeros(md.mesh.numberofelements,1); +%Basin separation TF +idb_tf = zeros(md.mesh.numberofelements,1); iid1 = find(md.mesh.x<=-1.6e6); for ii=1:md.mesh.numberofelements for vertex=1:3 if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1 - idbasin(ii) = 1; + idb_tf(ii) = 1; end end - if idbasin(ii)==0 %no vertex was found in basin 1 - idbasin(ii) = 2; + if idb_tf(ii)==0 %no vertex was found in basin 1 + idb_tf(ii) = 2; end end +% Basin separation default +idb_df = zeros(md.mesh.numberofelements,1); +iid1 = find(md.mesh.x<=-1.62e6); +for ii=1:md.mesh.numberofelements + for vertex=1:3 + if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1 + idb_df(ii) = 1; + end + end + if idb_df(ii)==0 %no vertex was found in basin 1 + idb_df(ii) = 2; + end +end +% Dimensionalities +nb_tf = 2; +nb_clv = 2; +nb_flmlt = 2; -%calving and frontalforcings parameters +%Calving parameters md.mask.ice_levelset = 1e4*(md.mask.ice_levelset + 0.5); -md.calving=calvingvonmises(); +md.calving.calvingrate = 0.3*ones(md.mesh.numberofvertices,1); +md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1); +md.levelset.migration_max = 10.0; %avoid fast advance/retreat of the front +%Frontal forcing parameters md.frontalforcings=frontalforcingsrignotautoregression(); -md.frontalforcings.num_basins = 2; -md.frontalforcings.basin_id = idbasin; +md.frontalforcings.num_basins = nb_tf; +md.frontalforcings.basin_id = idb_tf; md.frontalforcings.subglacial_discharge = 0.1*ones(md.mesh.numberofvertices,1); md.frontalforcings.beta0 = [0.05,0.01]; %intercept values of TF in basins [C] @@ -36,15 +56,26 @@ md.frontalforcings.phi = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters +% Floating Ice Melt parameters +md.basalforcings.floatingice_melting_rate = 0.1*ones(md.mesh.numberofvertices,1); + +% Covariance matrix +covtf = 1e-4*eye(nb_tf); +covclv = 1e-1*eye(nb_clv); +covclv(1,1) = 1/10*covclv(1,1); +covflmlt = 0.05*eye(nb_flmlt); +covglob = blkdiag(covtf,covclv,covflmlt); + %Stochastic forcing md.stochasticforcing.isstochasticforcing = 1; -md.stochasticforcing.fields = [{'FrontalForcingsRignotAutoregression'}]; -md.stochasticforcing.covariance = 1e-4*[[1.5,0.5];[0.5,0.4]]; %global covariance among- and between-fields +md.stochasticforcing.fields = [{'FrontalForcingsRignotAutoregression'},{'DefaultCalving'},{'FloatingMeltRate'}]; +md.stochasticforcing.defaultdimension = 2; +md.stochasticforcing.default_id = idb_df; +md.stochasticforcing.covariance = covglob; %global covariance among- and between-fields md.stochasticforcing.randomflag = 0; %determines true/false randomness -md.transient.ismovingfront = 1; -md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1); -md.levelset.migration_max = 1e10; +md.transient.ismovingfront = 1; +md.transient.isgroundingline = 1; -md.transient.requested_outputs = {'default','CalvingMeltingrate'}; +md.transient.requested_outputs = {'default','CalvingCalvingrate','CalvingMeltingrate','BasalforcingsFloatingiceMeltingRate'}; md.cluster=generic('name',oshostname(),'np',2); md=solve(md,'Transient'); @@ -52,12 +83,12 @@ %Fields and tolerances to track changes field_names ={... - 'Vx1' ,'Vy1' ,'Vel1' ,'Pressure1' ,'Bed1' ,'Surface1' ,'Thickness1' ,'MaskIceLevelset1', 'CalvingMeltingrate1',... - 'Vx2' ,'Vy2' ,'Vel2' ,'Pressure2' ,'Bed2' ,'Surface2' ,'Thickness2' ,'MaskIceLevelset2', 'CalvingMeltingrate2',... - 'Vx10','Vy10','Vel10','Pressure10','Bed10','Surface10','Thickness10','MaskIceLevelset10', 'CalvingMeltingrate10',... + 'Vx1' ,'Vy1' ,'Vel1' ,'Thickness1' ,'MaskIceLevelset1' ,'CalvingCalvingrate1' ,'CalvingMeltingrate1' ,'BasalforcingsFloatingiceMeltingRate1',... + 'Vx2' ,'Vy2' ,'Vel2' ,'Thickness2' ,'MaskIceLevelset2' ,'CalvingCalvingrate1' ,'CalvingMeltingrate2' ,'BasalforcingsFloatingiceMeltingRate2',... + 'Vx10','Vy10','Vel10','Thickness10','MaskIceLevelset10','CalvingCalvingrate10','CalvingMeltingrate10','BasalforcingsFloatingiceMeltingRate10',... }; field_tolerances={... - 1e-11,2e-11,2e-11,1e-12,2e-11,6e-12,9e-12,1e-12,1e-9,... - 2e-11,1e-11,1e-11,9e-12,2e-1,2e-11,2e-11,1e-11,1e-9,... - 2e-10,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10,1e-9,1e-9,... + 1e-11,2e-11,2e-11,1e-11,1e-9,1e-10,1e-10,1e-10,... + 2e-11,1e-11,1e-11,9e-11,2e-9,1e-10,1e-10,1e-10,... + 2e-10,1e-10,1e-10,1e-10,5e-9,1e-10,1e-10,1e-10,... }; field_values={... @@ -65,27 +96,24 @@ (md.results.TransientSolution(1).Vy),... (md.results.TransientSolution(1).Vel),... - (md.results.TransientSolution(1).Pressure),... - (md.results.TransientSolution(1).Base),... - (md.results.TransientSolution(1).Surface),... (md.results.TransientSolution(1).Thickness),... (md.results.TransientSolution(1).MaskIceLevelset),... + (md.results.TransientSolution(1).CalvingCalvingrate),... (md.results.TransientSolution(1).CalvingMeltingrate),... - (md.results.TransientSolution(2).Vx),... - (md.results.TransientSolution(2).Vy),... - (md.results.TransientSolution(2).Vel),... - (md.results.TransientSolution(2).Pressure),... - (md.results.TransientSolution(2).Base),... - (md.results.TransientSolution(2).Surface),... - (md.results.TransientSolution(2).Thickness),... - (md.results.TransientSolution(2).MaskIceLevelset),... - (md.results.TransientSolution(2).CalvingMeltingrate),... - (md.results.TransientSolution(10).Vx),... - (md.results.TransientSolution(10).Vy),... - (md.results.TransientSolution(10).Vel),... - (md.results.TransientSolution(10).Pressure),... - (md.results.TransientSolution(10).Base),... - (md.results.TransientSolution(10).Surface),... - (md.results.TransientSolution(10).Thickness),... - (md.results.TransientSolution(10).MaskIceLevelset),... - (md.results.TransientSolution(10).CalvingMeltingrate),... - }; + (md.results.TransientSolution(1).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(5).Vx),... + (md.results.TransientSolution(5).Vy),... + (md.results.TransientSolution(5).Vel),... + (md.results.TransientSolution(5).Thickness),... + (md.results.TransientSolution(5).MaskIceLevelset),... + (md.results.TransientSolution(5).CalvingCalvingrate),... + (md.results.TransientSolution(5).CalvingMeltingrate),... + (md.results.TransientSolution(5).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(10).Vx),... + (md.results.TransientSolution(10).Vy),... + (md.results.TransientSolution(10).Vel),... + (md.results.TransientSolution(10).Thickness),... + (md.results.TransientSolution(10).MaskIceLevelset),... + (md.results.TransientSolution(10).CalvingCalvingrate),... + (md.results.TransientSolution(10).CalvingMeltingrate),... + (md.results.TransientSolution(10).BasalforcingsFloatingiceMeltingRate),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test543.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.py (revision 26664) +++ /issm/trunk-jpl/test/NightlyRun/test543.py (revision 26665) @@ -1,6 +1,5 @@ -#Test Name: PigTranFrontalforcingsrignotautoregression +#Test Name: PigTranStochasticforcing import numpy as np -from calvingvonmises import * from frontalforcingsrignotautoregression import * from socket import gethostname @@ -19,23 +18,40 @@ md.timestepping.start_time = 0 md.timestepping.time_step = 1 -md.timestepping.final_time = 25 +md.timestepping.final_time = 10 -# Separate domain in 2 basins -idbasin = np.zeros((md.mesh.numberofelements,)) +# Basin separation TF +idb_tf = np.zeros((md.mesh.numberofelements,)) iid1 = np.where(md.mesh.x <= -1.6e6)[0] for ii in range(md.mesh.numberofelements): for vertex in range(3): if md.mesh.elements[ii][vertex] - 1 in iid1: # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing - idbasin[ii] = 1 + idb_tf[ii] = 1 if idbasin[ii] == 0: # no vertex was found in basin 1 for vertex in range(3): - idbasin[ii] = 2 + idb_tf[ii] = 2 +# Basin separation default +idb_df = np.zeros((md.mesh.numberofelements,)) +iid1 = np.where(md.mesh.x <= -1.62e6)[0] +for ii in range(md.mesh.numberofelements): + for vertex in range(3): + if md.mesh.elements[ii][vertex] - 1 in iid1: # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing + idb_df[ii] = 1 + if idbasin[ii] == 0: # no vertex was found in basin 1 + for vertex in range(3): + idb_df[ii] = 2 +#Dimensionalities +nb_tf = 2 +nb_clv = 2 +nb_flmlt = 2 -# Calving and frontalforcings parameters +# Calving parameters md.mask.ice_levelset = 1e4 * (md.mask.ice_levelset + 0.5) -md.calving = calvingvonmises() +md.calving.calvingrate = 0.3*np.ones((md.mesh.numberofvertices,)) +md.levelset.spclevelset = np.full((md.mesh.numberofvertices,), np.nan) +md.levelset.migration_max = 10.0 +#Frontal forcing parameters md.frontalforcings = frontalforcingsrignotautoregression() -md.frontalforcings.num_basins = 2 -md.frontalforcings.basin_id = idbasin - 1 # NOTE: offset because of 1-based vertex indexing +md.frontalforcings.num_basins = nb_tf +md.frontalforcings.basin_id = idb_tf - 1 # NOTE: offset because of 1-based vertex indexing md.frontalforcings.subglacial_discharge = 0.1 * np.ones((md.mesh.numberofvertices,)) md.frontalforcings.beta0 = np.array([[0.05, 0.01]]) # intercept values of TF in basins [C] @@ -45,18 +61,29 @@ md.frontalforcings.ar_timestep = 2 # timestep of the autoregressive model [yr] md.frontalforcings.phi = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]]) # autoregressive parameters +#Floating Ice Melt parameters +md.basalforcings.floatingice_melting_rate = 0.1*np.ones((md.mesh.numberofvertices,)) + +#Covariance matrix +covtf = 1e-4*np.identity(nb_tf) +covclv = 1e-1*np.identity(nb_clv) +covclv[0,0] = 1/10*covclv[0,0] +covflmlt = 0.05*np.identity(nb_flmlt) +covglob = np.zeros([6,6]) +covglob[0:2,0:2] = covtf +covglob[2:4,2:4] = covclv +covglob[4:6,4:6] = covflmlt # Stochastic forcing md.stochasticforcing.isstochasticforcing = 1 -md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression'] +md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression','DefaultCalving','FloatingMeltRate'] +md.stochasticforcing.defauldimension = 2 +md.stochasticforcing.default_id = idb_df +md.stochasticforcing.covariance = covglob # global covariance among- and between-fields +md.stochasticforcing.randomflag = 0 # determines true/false randomness -md.stochasticforcing.covariance = 1e-4 * np.array([[1.5, 0.5], [0.5, 0.4]]) # global covariance among- and between-fields -md.stochasticforcing.randomflag = 0 # determines true/false randomness +md.transient.ismovingfront = 1 +md.transient.isgroundingline = 1 - -md.transient.ismovingfront = 1 -md.levelset.spclevelset = np.full((md.mesh.numberofvertices,), np.nan) -md.levelset.migration_max = 1e10 - -md.transient.requested_outputs = ['default', 'CalvingMeltingrate'] +md.transient.requested_outputs = ['default', 'CalvingCalvingrate','CalvingMeltingrate','BasalforcingsFloatingiceMeltingRate'] md.cluster = generic('name', gethostname(), 'np', 2) md = solve(md, 'Transient') @@ -64,12 +91,13 @@ # Fields and tolerances to track changes field_names = [ - 'Vx1', 'Vy1', 'Vel1', 'Pressure1', 'Bed1', 'Surface1', 'Thickness1', 'MaskIceLevelset1', 'CalvingMetlingRate1', - 'Vx2', 'Vy2', 'Vel2', 'Pressure2', 'Bed2', 'Surface2', 'Thickness2', 'MaskIceLevelset2', 'CalvingMetlingRate2', - 'Vx10', 'Vy10', 'Vel10', 'Pressure10', 'Bed10', 'Surface10', 'Thickness10', 'MaskIceLevelset10', 'CalvingMetlingRate10' + 'Vx1' ,'Vy1' ,'Vel1' ,'Thickness1' ,'MaskIceLevelset1' ,'CalvingCalvingrate1' ,'CalvingMeltingrate1' ,'BasalforcingsFloatingiceMeltingRate1', + 'Vx2' ,'Vy2' ,'Vel2' ,'Thickness2' ,'MaskIceLevelset2' ,'CalvingCalvingrate1' ,'CalvingMeltingrate2' ,'BasalforcingsFloatingiceMeltingRate2', + 'Vx10','Vy10','Vel10','Thickness10','MaskIceLevelset10','CalvingCalvingrate10','CalvingMeltingrate10','BasalforcingsFloatingiceMeltingRate10' ] + field_tolerances = [ - 1e-11, 2e-11, 2e-11, 1e-12, 2e-11, 6e-12, 9e-12, 1e-12, 1e-9, - 2e-11, 1e-11, 1e-11, 9e-12, 2e-1, 2e-11, 2e-11, 1e-11, 1e-9, - 2e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-9, 1e-9 + 1e-11,2e-11,2e-11,1e-11,1e-9,1e-10,1e-10,1e-10, + 2e-11,1e-11,1e-11,9e-11,2e-9,1e-10,1e-10,1e-10, + 2e-10,1e-10,1e-10,1e-10,5e-9,1e-10,1e-10,1e-10 ] field_values = [ @@ -77,27 +105,24 @@ md.results.TransientSolution[0].Vy, md.results.TransientSolution[0].Vel, - md.results.TransientSolution[0].Pressure, - md.results.TransientSolution[0].Base, - md.results.TransientSolution[0].Surface, md.results.TransientSolution[0].Thickness, md.results.TransientSolution[0].MaskIceLevelset, + md.results.TransientSolution[0].CalvingCalvingrate, md.results.TransientSolution[0].CalvingMeltingrate, - md.results.TransientSolution[1].Vx, - md.results.TransientSolution[1].Vy, - md.results.TransientSolution[1].Vel, - md.results.TransientSolution[1].Pressure, - md.results.TransientSolution[1].Base, - md.results.TransientSolution[1].Surface, - md.results.TransientSolution[1].Thickness, - md.results.TransientSolution[1].MaskIceLevelset, - md.results.TransientSolution[1].CalvingMeltingrate, + md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[4].Vx, + md.results.TransientSolution[4].Vy, + md.results.TransientSolution[4].Vel, + md.results.TransientSolution[4].Thickness, + md.results.TransientSolution[4].MaskIceLevelset, + md.results.TransientSolution[4].CalvingCalvingrate, + md.results.TransientSolution[4].CalvingMeltingrate, + md.results.TransientSolution[4].BasalforcingsFloatingiceMeltingRate, md.results.TransientSolution[9].Vx, md.results.TransientSolution[9].Vy, md.results.TransientSolution[9].Vel, - md.results.TransientSolution[9].Pressure, - md.results.TransientSolution[9].Base, - md.results.TransientSolution[9].Surface, md.results.TransientSolution[9].Thickness, md.results.TransientSolution[9].MaskIceLevelset, - md.results.TransientSolution[9].CalvingMeltingrate + md.results.TransientSolution[9].CalvingCalvingrate, + md.results.TransientSolution[9].CalvingMeltingrate, + md.results.TransientSolution[9].BasalforcingsFloatingiceMeltingRate ]