source: issm/oecreview/Archive/25834-26739/ISSM-26664-26665.diff

Last change on this file was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739

File size: 15.7 KB
  • ../trunk-jpl/test/NightlyRun/test543.m

    Cannot display: file marked as a binary type.
    svn:mime-type = application/octet-stream
     
    1 %Test Name: PigTranFrontalforcingsrignotautoregression
     1%Test Name: PigTranStochasticforcing
    22md=triangle(model(),'../Exp/Pig.exp',10000.);
    33md=setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp');
    44md=parameterize(md,'../Par/Pig.par');
     
    55md=setflowequation(md,'SSA','all');
    66md.timestepping.start_time = 0;
    77md.timestepping.time_step  = 1;
    8 md.timestepping.final_time = 25;
     8md.timestepping.final_time = 10;
    99
    10 %separate domain in 2 basins
    11 idbasin = zeros(md.mesh.numberofelements,1);
     10%Basin separation TF
     11idb_tf = zeros(md.mesh.numberofelements,1);
    1212iid1    = find(md.mesh.x<=-1.6e6);
    1313for ii=1:md.mesh.numberofelements
    1414    for vertex=1:3
    1515        if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1
    16             idbasin(ii) = 1;
     16            idb_tf(ii) = 1;
    1717        end
    1818    end
    19     if idbasin(ii)==0 %no vertex was found in basin 1
    20         idbasin(ii) = 2;
     19    if idb_tf(ii)==0 %no vertex was found in basin 1
     20        idb_tf(ii) = 2;
    2121    end
    2222end
     23% Basin separation default
     24idb_df = zeros(md.mesh.numberofelements,1);
     25iid1   = find(md.mesh.x<=-1.62e6);
     26for ii=1:md.mesh.numberofelements
     27    for vertex=1:3
     28        if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1
     29            idb_df(ii) = 1;
     30        end
     31    end
     32    if idb_df(ii)==0 %no vertex was found in basin 1
     33        idb_df(ii) = 2;
     34    end
     35end
     36% Dimensionalities
     37nb_tf    = 2;
     38nb_clv   = 2;
     39nb_flmlt = 2;
    2340
    24 %calving and frontalforcings parameters
     41%Calving parameters
    2542md.mask.ice_levelset = 1e4*(md.mask.ice_levelset + 0.5);
    26 md.calving=calvingvonmises();
     43md.calving.calvingrate = 0.3*ones(md.mesh.numberofvertices,1);
     44md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1);
     45md.levelset.migration_max = 10.0; %avoid fast advance/retreat of the front
     46%Frontal forcing parameters
    2747md.frontalforcings=frontalforcingsrignotautoregression();
    28 md.frontalforcings.num_basins = 2;
    29 md.frontalforcings.basin_id = idbasin;
     48md.frontalforcings.num_basins           = nb_tf;
     49md.frontalforcings.basin_id             = idb_tf;
    3050md.frontalforcings.subglacial_discharge = 0.1*ones(md.mesh.numberofvertices,1);
    3151md.frontalforcings.beta0                = [0.05,0.01]; %intercept values of TF in basins [C]
    3252md.frontalforcings.beta1                = [0.001,0.0001]; %trend values of TF in basins [C/yr]
     
    3555md.frontalforcings.ar_timestep          = 2; %timestep of the autoregressive model [yr]
    3656md.frontalforcings.phi                  = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters
    3757
     58% Floating Ice Melt parameters
     59md.basalforcings.floatingice_melting_rate = 0.1*ones(md.mesh.numberofvertices,1);
     60
     61% Covariance matrix
     62covtf       = 1e-4*eye(nb_tf);
     63covclv      = 1e-1*eye(nb_clv);
     64covclv(1,1) = 1/10*covclv(1,1);
     65covflmlt    = 0.05*eye(nb_flmlt);
     66covglob     = blkdiag(covtf,covclv,covflmlt);
     67
    3868%Stochastic forcing
    3969md.stochasticforcing.isstochasticforcing = 1;
    40 md.stochasticforcing.fields              = [{'FrontalForcingsRignotAutoregression'}];
    41 md.stochasticforcing.covariance          = 1e-4*[[1.5,0.5];[0.5,0.4]]; %global covariance among- and between-fields
     70md.stochasticforcing.fields              = [{'FrontalForcingsRignotAutoregression'},{'DefaultCalving'},{'FloatingMeltRate'}];
     71md.stochasticforcing.defaultdimension    = 2;
     72md.stochasticforcing.default_id          = idb_df;
     73md.stochasticforcing.covariance          = covglob; %global covariance among- and between-fields
    4274md.stochasticforcing.randomflag          = 0; %determines true/false randomness
    4375
    44 md.transient.ismovingfront = 1;
    45 md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1);
    46 md.levelset.migration_max = 1e10;
     76md.transient.ismovingfront   = 1;
     77md.transient.isgroundingline = 1;
    4778
    48 md.transient.requested_outputs = {'default','CalvingMeltingrate'};
     79md.transient.requested_outputs = {'default','CalvingCalvingrate','CalvingMeltingrate','BasalforcingsFloatingiceMeltingRate'};
    4980md.cluster=generic('name',oshostname(),'np',2);
    5081md=solve(md,'Transient');
    5182
    5283%Fields and tolerances to track changes
    5384field_names ={...
    54    'Vx1' ,'Vy1' ,'Vel1' ,'Pressure1' ,'Bed1' ,'Surface1' ,'Thickness1' ,'MaskIceLevelset1', 'CalvingMeltingrate1',...
    55    'Vx2' ,'Vy2' ,'Vel2' ,'Pressure2' ,'Bed2' ,'Surface2' ,'Thickness2' ,'MaskIceLevelset2', 'CalvingMeltingrate2',...
    56    'Vx10','Vy10','Vel10','Pressure10','Bed10','Surface10','Thickness10','MaskIceLevelset10', 'CalvingMeltingrate10',...
     85   'Vx1' ,'Vy1' ,'Vel1' ,'Thickness1' ,'MaskIceLevelset1' ,'CalvingCalvingrate1' ,'CalvingMeltingrate1' ,'BasalforcingsFloatingiceMeltingRate1',...
     86   'Vx2' ,'Vy2' ,'Vel2' ,'Thickness2' ,'MaskIceLevelset2' ,'CalvingCalvingrate1' ,'CalvingMeltingrate2' ,'BasalforcingsFloatingiceMeltingRate2',...
     87   'Vx10','Vy10','Vel10','Thickness10','MaskIceLevelset10','CalvingCalvingrate10','CalvingMeltingrate10','BasalforcingsFloatingiceMeltingRate10',...
    5788   };
    5889field_tolerances={...
    59    1e-11,2e-11,2e-11,1e-12,2e-11,6e-12,9e-12,1e-12,1e-9,...
    60    2e-11,1e-11,1e-11,9e-12,2e-1,2e-11,2e-11,1e-11,1e-9,...
    61    2e-10,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10,1e-9,1e-9,...
     90   1e-11,2e-11,2e-11,1e-11,1e-9,1e-10,1e-10,1e-10,...
     91   2e-11,1e-11,1e-11,9e-11,2e-9,1e-10,1e-10,1e-10,...
     92   2e-10,1e-10,1e-10,1e-10,5e-9,1e-10,1e-10,1e-10,...
    6293   };
    6394field_values={...
    6495   (md.results.TransientSolution(1).Vx),...
    6596   (md.results.TransientSolution(1).Vy),...
    6697   (md.results.TransientSolution(1).Vel),...
    67    (md.results.TransientSolution(1).Pressure),...
    68    (md.results.TransientSolution(1).Base),...
    69    (md.results.TransientSolution(1).Surface),...
    7098   (md.results.TransientSolution(1).Thickness),...
    7199   (md.results.TransientSolution(1).MaskIceLevelset),...
     100   (md.results.TransientSolution(1).CalvingCalvingrate),...
    72101   (md.results.TransientSolution(1).CalvingMeltingrate),...
    73    (md.results.TransientSolution(2).Vx),...
    74    (md.results.TransientSolution(2).Vy),...
    75    (md.results.TransientSolution(2).Vel),...
    76    (md.results.TransientSolution(2).Pressure),...
    77    (md.results.TransientSolution(2).Base),...
    78    (md.results.TransientSolution(2).Surface),...
    79    (md.results.TransientSolution(2).Thickness),...
    80    (md.results.TransientSolution(2).MaskIceLevelset),...
    81    (md.results.TransientSolution(2).CalvingMeltingrate),...
    82    (md.results.TransientSolution(10).Vx),...
    83    (md.results.TransientSolution(10).Vy),...
    84    (md.results.TransientSolution(10).Vel),...
    85    (md.results.TransientSolution(10).Pressure),...
    86    (md.results.TransientSolution(10).Base),...
    87    (md.results.TransientSolution(10).Surface),...
    88    (md.results.TransientSolution(10).Thickness),...
    89    (md.results.TransientSolution(10).MaskIceLevelset),...
    90    (md.results.TransientSolution(10).CalvingMeltingrate),...
    91    };
     102   (md.results.TransientSolution(1).BasalforcingsFloatingiceMeltingRate),...
     103   (md.results.TransientSolution(5).Vx),...
     104   (md.results.TransientSolution(5).Vy),...
     105   (md.results.TransientSolution(5).Vel),...
     106   (md.results.TransientSolution(5).Thickness),...
     107   (md.results.TransientSolution(5).MaskIceLevelset),...
     108   (md.results.TransientSolution(5).CalvingCalvingrate),...
     109   (md.results.TransientSolution(5).CalvingMeltingrate),...
     110   (md.results.TransientSolution(5).BasalforcingsFloatingiceMeltingRate),...
     111        (md.results.TransientSolution(10).Vx),...
     112        (md.results.TransientSolution(10).Vy),...
     113        (md.results.TransientSolution(10).Vel),...
     114        (md.results.TransientSolution(10).Thickness),...
     115        (md.results.TransientSolution(10).MaskIceLevelset),...
     116        (md.results.TransientSolution(10).CalvingCalvingrate),...
     117        (md.results.TransientSolution(10).CalvingMeltingrate),...
     118        (md.results.TransientSolution(10).BasalforcingsFloatingiceMeltingRate),...
     119        };
  • ../trunk-jpl/test/NightlyRun/test543.py

     
    1 #Test Name: PigTranFrontalforcingsrignotautoregression
     1#Test Name: PigTranStochasticforcing
    22import numpy as np
    33
    4 from calvingvonmises import *
    54from frontalforcingsrignotautoregression import *
    65from socket import gethostname
    76from model import *
     
    1817md = setflowequation(md, 'SSA', 'all')
    1918md.timestepping.start_time = 0
    2019md.timestepping.time_step = 1
    21 md.timestepping.final_time = 25
     20md.timestepping.final_time = 10
    2221
    23 # Separate domain in 2 basins
    24 idbasin = np.zeros((md.mesh.numberofelements,))
     22# Basin separation TF
     23idb_tf = np.zeros((md.mesh.numberofelements,))
    2524iid1 = np.where(md.mesh.x <= -1.6e6)[0]
    2625for ii in range(md.mesh.numberofelements):
    2726    for vertex in range(3):
    2827        if md.mesh.elements[ii][vertex] - 1 in iid1:  # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing
    29             idbasin[ii] = 1
     28            idb_tf[ii] = 1
    3029    if idbasin[ii] == 0:  # no vertex was found in basin 1
    3130        for vertex in range(3):
    32             idbasin[ii] = 2
     31            idb_tf[ii] = 2
     32# Basin separation default
     33idb_df = np.zeros((md.mesh.numberofelements,))
     34iid1 = np.where(md.mesh.x <= -1.62e6)[0]
     35for ii in range(md.mesh.numberofelements):
     36    for vertex in range(3):
     37        if md.mesh.elements[ii][vertex] - 1 in iid1:  # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing
     38            idb_df[ii] = 1
     39    if idbasin[ii] == 0:  # no vertex was found in basin 1
     40        for vertex in range(3):
     41            idb_df[ii] = 2
     42#Dimensionalities
     43nb_tf    = 2
     44nb_clv   = 2
     45nb_flmlt = 2
    3346
    34 # Calving and frontalforcings parameters
     47# Calving parameters
    3548md.mask.ice_levelset = 1e4 * (md.mask.ice_levelset + 0.5)
    36 md.calving = calvingvonmises()
     49md.calving.calvingrate = 0.3*np.ones((md.mesh.numberofvertices,))
     50md.levelset.spclevelset = np.full((md.mesh.numberofvertices,), np.nan)
     51md.levelset.migration_max = 10.0
     52#Frontal forcing parameters
    3753md.frontalforcings = frontalforcingsrignotautoregression()
    38 md.frontalforcings.num_basins = 2
    39 md.frontalforcings.basin_id = idbasin - 1  # NOTE: offset because of 1-based vertex indexing
     54md.frontalforcings.num_basins = nb_tf
     55md.frontalforcings.basin_id = idb_tf - 1  # NOTE: offset because of 1-based vertex indexing
    4056md.frontalforcings.subglacial_discharge = 0.1 * np.ones((md.mesh.numberofvertices,))
    4157md.frontalforcings.beta0 = np.array([[0.05, 0.01]])  # intercept values of TF in basins [C]
    4258md.frontalforcings.beta1 = np.array([[0.001, 0.0001]])  # trend values of TF in basins [C/yr]
     
    4460md.frontalforcings.ar_order = 4
    4561md.frontalforcings.ar_timestep = 2  # timestep of the autoregressive model [yr]
    4662md.frontalforcings.phi = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]])  # autoregressive parameters
     63#Floating Ice Melt parameters
     64md.basalforcings.floatingice_melting_rate = 0.1*np.ones((md.mesh.numberofvertices,))
    4765
     66#Covariance matrix
     67covtf            = 1e-4*np.identity(nb_tf)
     68covclv           = 1e-1*np.identity(nb_clv)
     69covclv[0,0]      = 1/10*covclv[0,0]
     70covflmlt         = 0.05*np.identity(nb_flmlt)
     71covglob          = np.zeros([6,6])
     72covglob[0:2,0:2] = covtf
     73covglob[2:4,2:4] = covclv
     74covglob[4:6,4:6] = covflmlt
     75
    4876# Stochastic forcing
    4977md.stochasticforcing.isstochasticforcing = 1
    50 md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression']
     78md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression','DefaultCalving','FloatingMeltRate']
     79md.stochasticforcing.defauldimension = 2
     80md.stochasticforcing.default_id      = idb_df
     81md.stochasticforcing.covariance      = covglob # global covariance among- and between-fields
     82md.stochasticforcing.randomflag      = 0 # determines true/false randomness
    5183
    52 md.stochasticforcing.covariance = 1e-4 * np.array([[1.5, 0.5], [0.5, 0.4]]) # global covariance among- and between-fields
    53 md.stochasticforcing.randomflag = 0 # determines true/false randomness
     84md.transient.ismovingfront   = 1
     85md.transient.isgroundingline = 1
    5486
    55 
    56 md.transient.ismovingfront = 1
    57 md.levelset.spclevelset = np.full((md.mesh.numberofvertices,), np.nan)
    58 md.levelset.migration_max = 1e10
    59 
    60 md.transient.requested_outputs = ['default', 'CalvingMeltingrate']
     87md.transient.requested_outputs = ['default', 'CalvingCalvingrate','CalvingMeltingrate','BasalforcingsFloatingiceMeltingRate']
    6188md.cluster = generic('name', gethostname(), 'np', 2)
    6289md = solve(md, 'Transient')
    6390
    6491# Fields and tolerances to track changes
    6592field_names = [
    66     'Vx1', 'Vy1', 'Vel1', 'Pressure1', 'Bed1', 'Surface1', 'Thickness1', 'MaskIceLevelset1', 'CalvingMetlingRate1',
    67     'Vx2', 'Vy2', 'Vel2', 'Pressure2', 'Bed2', 'Surface2', 'Thickness2', 'MaskIceLevelset2', 'CalvingMetlingRate2',
    68     'Vx10', 'Vy10', 'Vel10', 'Pressure10', 'Bed10', 'Surface10', 'Thickness10', 'MaskIceLevelset10', 'CalvingMetlingRate10'
     93    'Vx1' ,'Vy1' ,'Vel1' ,'Thickness1' ,'MaskIceLevelset1' ,'CalvingCalvingrate1' ,'CalvingMeltingrate1' ,'BasalforcingsFloatingiceMeltingRate1',
     94    'Vx2' ,'Vy2' ,'Vel2' ,'Thickness2' ,'MaskIceLevelset2' ,'CalvingCalvingrate1' ,'CalvingMeltingrate2' ,'BasalforcingsFloatingiceMeltingRate2',
     95    'Vx10','Vy10','Vel10','Thickness10','MaskIceLevelset10','CalvingCalvingrate10','CalvingMeltingrate10','BasalforcingsFloatingiceMeltingRate10'
    6996]
     97
    7098field_tolerances = [
    71     1e-11, 2e-11, 2e-11, 1e-12, 2e-11, 6e-12, 9e-12, 1e-12, 1e-9,
    72     2e-11, 1e-11, 1e-11, 9e-12, 2e-1, 2e-11, 2e-11, 1e-11, 1e-9,
    73     2e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-9, 1e-9
     99    1e-11,2e-11,2e-11,1e-11,1e-9,1e-10,1e-10,1e-10,
     100    2e-11,1e-11,1e-11,9e-11,2e-9,1e-10,1e-10,1e-10,
     101    2e-10,1e-10,1e-10,1e-10,5e-9,1e-10,1e-10,1e-10
    74102]
    75103field_values = [
    76104    md.results.TransientSolution[0].Vx,
    77105    md.results.TransientSolution[0].Vy,
    78106    md.results.TransientSolution[0].Vel,
    79     md.results.TransientSolution[0].Pressure,
    80     md.results.TransientSolution[0].Base,
    81     md.results.TransientSolution[0].Surface,
    82107    md.results.TransientSolution[0].Thickness,
    83108    md.results.TransientSolution[0].MaskIceLevelset,
     109    md.results.TransientSolution[0].CalvingCalvingrate,
    84110    md.results.TransientSolution[0].CalvingMeltingrate,
    85     md.results.TransientSolution[1].Vx,
    86     md.results.TransientSolution[1].Vy,
    87     md.results.TransientSolution[1].Vel,
    88     md.results.TransientSolution[1].Pressure,
    89     md.results.TransientSolution[1].Base,
    90     md.results.TransientSolution[1].Surface,
    91     md.results.TransientSolution[1].Thickness,
    92     md.results.TransientSolution[1].MaskIceLevelset,
    93     md.results.TransientSolution[1].CalvingMeltingrate,
     111    md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate,
     112    md.results.TransientSolution[4].Vx,
     113    md.results.TransientSolution[4].Vy,
     114    md.results.TransientSolution[4].Vel,
     115    md.results.TransientSolution[4].Thickness,
     116    md.results.TransientSolution[4].MaskIceLevelset,
     117    md.results.TransientSolution[4].CalvingCalvingrate,
     118    md.results.TransientSolution[4].CalvingMeltingrate,
     119    md.results.TransientSolution[4].BasalforcingsFloatingiceMeltingRate,
    94120    md.results.TransientSolution[9].Vx,
    95121    md.results.TransientSolution[9].Vy,
    96122    md.results.TransientSolution[9].Vel,
    97     md.results.TransientSolution[9].Pressure,
    98     md.results.TransientSolution[9].Base,
    99     md.results.TransientSolution[9].Surface,
    100123    md.results.TransientSolution[9].Thickness,
    101124    md.results.TransientSolution[9].MaskIceLevelset,
    102     md.results.TransientSolution[9].CalvingMeltingrate
     125    md.results.TransientSolution[9].CalvingCalvingrate,
     126    md.results.TransientSolution[9].CalvingMeltingrate,
     127    md.results.TransientSolution[9].BasalforcingsFloatingiceMeltingRate
    103128]
Note: See TracBrowser for help on using the repository browser.