Changeset 909


Ignore:
Timestamp:
06/11/09 14:56:29 (15 years ago)
Author:
seroussi
Message:

simplified parameter file

Location:
issm/trunk/test/Validation/MacAyealVsPattyn
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/test/Validation/MacAyealVsPattyn/test1_iceshelf/Square.par

    r16 r909  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.8;
    6         md.rho_ice=917;
    7         md.rho_water=1023;
    8         di=md.rho_ice/md.rho_water;
    9         md.yts=365*24*3600;
    10         md.heatcapacity=2009;
    11         md.thermalconductivity=2.2; %W/mK
    12         md.beta=9.8*10^-8;
     3disp('      creating thickness');
     4hmin=300;
     5hmax=1000;
     6ymin=min(md.y);
     7ymax=max(md.y);
     8md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
     9md.firn_layer=10*ones(md.numberofgrids,1);
     10md.bed=-md.rho_ice/md.rho_water*md.thickness;
     11md.surface=md.bed+md.thickness;
    1312
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     13disp('      creating temperature');
     14md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    2015
    21         %statics
    22         md.eps_rel=0.01;
    23         md.eps_abs=10;
    24         md.lowmem=1;
    25         if md.numberofgrids<1000000,
    26         md.sparsity=.001;
    27         else
    28         md.sparsity=.0001;
    29         end
     16disp('      creating drag');
     17md.drag_type=2; %0 none 1 plastic 2 viscous
     18md.drag=200*ones(md.numberofgrids,1); %q=1.
     19%Take care of iceshelves: no basal drag
     20pos=find(md.elementoniceshelf);
     21md.drag(md.elements(pos,:))=0;
     22md.p=ones(md.numberofelements,1);
     23md.q=ones(md.numberofelements,1);
    3024
    31         %dynamics
    32         md.dt=1*md.yts; %1 year
    33         md.ndt=md.dt*10;
    34         md.artificial_diffusivity=1;
     25disp('      creating flow law paramter');
     26md.B=paterson(md.observed_temperature);
     27md.n=3*ones(md.numberofelements,1);
    3528
    36         %control
    37         md.control_type={'drag'}; %'drag', 'B'
    38         md.nsteps=5;
    39         md.tolx=10^-4;
    40         md.maxiter=20;
    41         md.optscal=10;
    42         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    43         md.meanvel=1000/md.yts; %1000 meters/year
    44         md.epsvel=eps;
     29disp('      boundary conditions for diagnostic model: ');
     30md=SetMarineIceSheetBC(md,'Front.exp');
    4531
    46 
    47         disp('      creating thickness');
    48         hmin=300;
    49         hmax=1000;
    50         ymin=min(md.y);
    51         ymax=max(md.y);
    52         md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
    53         md.firn_layer=10*ones(md.numberofgrids,1);
    54         md.bed=-di*md.thickness+10;
    55         md.surface=md.bed+md.thickness;
    56 
    57         disp('      creating velocities');
    58         md.vx_obs=zeros(md.numberofgrids,1);
    59         md.vy_obs=zeros(md.numberofgrids,1);
    60         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    61        
    62         disp('      creating drag');
    63         md.drag_type=2; %0 none 1 plastic 2 viscous
    64         md.drag=200*ones(md.numberofgrids,1); %q=1.
    65         %Take care of iceshelves: no basal drag
    66         pos=find(md.elementoniceshelf);
    67         md.drag(md.elements(pos,:))=0;
    68         md.p=ones(md.numberofelements,1);
    69         md.q=ones(md.numberofelements,1);
    70 
    71         disp('      creating temperature');
    72         md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    73 
    74         disp('      creating flow law paramter');
    75         md.B=paterson(md.observed_temperature);
    76         md.n=3*ones(md.numberofelements,1);
    77 
    78         disp('      creating accumulation rates');
    79         md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
    80         md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
    81 
    82         %Deal with boundary conditions:
    83        
    84         disp('      boundary conditions for diagnostic model: ');
    85         %Build gridonicefront, array of boundary grids belonging to the icefront:
    86         gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    87         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    88 
    89         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    90         pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
    91         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    92 
    93         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    94         md.segmentonneumann_diag=md.segments(pos,:);
    95         md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    96 
    97         disp('      boundary conditions for prognostic model: ');
    98         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    99         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    100         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    101         md.segmentonneumann_prog=md.segments(pos,:);
    102         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    103         md.neumannvalues_prog(:)=NaN; %free radiation
    104        
    105         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    106         md.segmentonneumann_prog2=md.segments(pos,:);
    107         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    108         md.neumannvalues_prog2(:)=NaN; %free radiation
    109        
    110         disp('      boundary conditions for thermal model: ');
    111         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    112         md.dirichletvalues_thermal=md.observed_temperature;
    113         md.geothermalflux=zeros(md.numberofgrids,1);
    114         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    115 
    116 
    117        
    118 
    119 % Some Cielo code, ignore.
    120 if strcmp(md.cluster,'yes')
    121         ServerDisconnect;
    122 end   
     32%Parallel options
     33md.np=3;
     34md.time=50;
     35md.waitonlock=1;

  • issm/trunk/test/Validation/MacAyealVsPattyn/test1_iceshelf/runme.m

    r899 r909  
    2424for i=1:md.numberofgrids2d
    2525        for j=1:(md.numlayers-1)
    26                 grid_vel=grid_vel+1/(2*(md.numlayers-1))*(md.results.diagnostic.veli+j*md.numberofgrids2d,1)+md.results.diagnostic.vel(i+(j-1)*md.numberofgrids2d,1));
     26                grid_vel=grid_vel+1/(2*(md.numlayers-1))*(md.results.diagnostic.vel(i+j*md.numberofgrids2d,1)+md.results.diagnostic.vel(i+(j-1)*md.numberofgrids2d,1));
    2727        end
    2828        vel_3d(i,1)=grid_vel;

  • issm/trunk/test/Validation/MacAyealVsPattyn/test2_icesheet/Square.par

    r16 r909  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.8;
    6         md.rho_ice=917;
    7         md.rho_water=1023;
    8         di=md.rho_ice/md.rho_water;
    9         md.yts=365*24*3600;
    10         md.heatcapacity=2009;
    11         md.thermalconductivity=2.2; %W/mK
    12         md.beta=9.8*10^-8;
     3disp('      creating thickness');
     4hmin=300;
     5hmax=1000;
     6ymin=min(md.y);
     7ymax=max(md.y);
     8md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
     9md.firn_layer=10*ones(md.numberofgrids,1);
     10md.bed=-md.rho_ice/md.rho_water*md.thickness+10;
     11md.surface=md.bed+md.thickness;
    1312
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     13disp('      creating temperature');
     14md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    2015
    21         %statics
    22         md.eps_rel=0.01;
    23         md.eps_abs=10;
    24         md.lowmem=1;
    25         if md.numberofgrids<1000000,
    26         md.sparsity=.001;
    27         else
    28         md.sparsity=.0001;
    29         end
     16disp('      creating drag');
     17md.drag_type=2; %0 none 1 plastic 2 viscous
     18md.drag=10*ones(md.numberofgrids,1); %q=1.
     19%Take care of iceshelves: no basal drag
     20pos=find(md.elementoniceshelf);
     21md.drag(md.elements(pos,:))=0;
     22md.p=ones(md.numberofelements,1);
     23md.q=ones(md.numberofelements,1);
    3024
    31         %dynamics
    32         md.dt=1*md.yts; %1 year
    33         md.ndt=md.dt*10;
    34         md.artificial_diffusivity=1;
     25disp('      creating flow law paramter');
     26md.B=paterson(md.observed_temperature);
     27md.n=3*ones(md.numberofelements,1);
    3528
    36         %control
    37         md.control_type={'drag'}; %'drag', 'B'
    38         md.nsteps=5;
    39         md.tolx=10^-4;
    40         md.maxiter=20;
    41         md.optscal=10;
    42         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    43         md.meanvel=1000/md.yts; %1000 meters/year
    44         md.epsvel=eps;
     29%Deal with boundary conditions:
     30disp('      boundary conditions for diagnostic model: ');
     31md=SetMarineIceSheetBC(md,'Front.exp');
    4532
    46 
    47         disp('      creating thickness');
    48         hmin=300;
    49         hmax=1000;
    50         ymin=min(md.y);
    51         ymax=max(md.y);
    52         md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
    53         md.firn_layer=10*ones(md.numberofgrids,1);
    54         md.bed=-di*md.thickness+10;
    55         md.surface=md.bed+md.thickness;
    56 
    57         disp('      creating velocities');
    58         md.vx_obs=zeros(md.numberofgrids,1);
    59         md.vy_obs=zeros(md.numberofgrids,1);
    60         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    61        
    62         disp('      creating drag');
    63         md.drag_type=2; %0 none 1 plastic 2 viscous
    64         md.drag=10*ones(md.numberofgrids,1); %q=1.
    65         %Take care of iceshelves: no basal drag
    66         pos=find(md.elementoniceshelf);
    67         md.drag(md.elements(pos,:))=0;
    68         md.p=ones(md.numberofelements,1);
    69         md.q=ones(md.numberofelements,1);
    70 
    71         disp('      creating temperature');
    72         md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    73 
    74         disp('      creating flow law paramter');
    75         md.B=paterson(md.observed_temperature);
    76         md.n=3*ones(md.numberofelements,1);
    77 
    78         disp('      creating accumulation rates');
    79         md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
    80         md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
    81 
    82         %Deal with boundary conditions:
    83        
    84         disp('      boundary conditions for diagnostic model: ');
    85         %Build gridonicefront, array of boundary grids belonging to the icefront:
    86         gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    87         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    88 
    89         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    90         pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
    91         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    92 
    93         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    94         md.segmentonneumann_diag=md.segments(pos,:);
    95         md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    96 
    97         disp('      boundary conditions for prognostic model: ');
    98         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    99         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    100         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    101         md.segmentonneumann_prog=md.segments(pos,:);
    102         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    103         md.neumannvalues_prog(:)=NaN; %free radiation
    104        
    105         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    106         md.segmentonneumann_prog2=md.segments(pos,:);
    107         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    108         md.neumannvalues_prog2(:)=NaN; %free radiation
    109        
    110         disp('      boundary conditions for thermal model: ');
    111         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    112         md.dirichletvalues_thermal=md.observed_temperature;
    113         md.geothermalflux=zeros(md.numberofgrids,1);
    114         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    115 
    116 
    117        
    118 
    119 % Some Cielo code, ignore.
    120 if strcmp(md.cluster,'yes')
    121         ServerDisconnect;
    122 end   
     33%Parallel options
     34md.np=3;
     35md.time=50;
     36md.waitonlock=1;

  • issm/trunk/test/Validation/MacAyealVsPattyn/test3_icesheet_iceshelf/Square.par

    r16 r909  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.8;
    6         md.rho_ice=917;
    7         md.rho_water=1023;
    8         di=md.rho_ice/md.rho_water;
    9         md.yts=365*24*3600;
    10         md.heatcapacity=2009;
    11         md.thermalconductivity=2.2; %W/mK
    12         md.beta=9.8*10^-8;
     3disp('      creating thickness'); hmin=300;
     4hmax=1000;
     5ymin=min(md.y);
     6ymax=max(md.y);
     7md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
     8md.firn_layer=10*ones(md.numberofgrids,1);
     9md.bed=-md.rho_ice/md.rho_water*md.thickness+10;
     10md.surface=md.bed+md.thickness;
    1311
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     12disp('      creating drag');
     13md.drag_type=2; %0 none 1 plastic 2 viscous
     14md.drag=200*ones(md.numberofgrids,1); %q=1.
     15%Take care of iceshelves: no basal drag
     16pos=find(md.elementoniceshelf);
     17md.drag(md.elements(pos,:))=0;
     18md.p=ones(md.numberofelements,1);
     19md.q=ones(md.numberofelements,1);
    2020
    21         %statics
    22         md.eps_rel=0.01;
    23         md.eps_abs=10;
    24         md.lowmem=1;
    25         if md.numberofgrids<1000000,
    26         md.sparsity=.001;
    27         else
    28         md.sparsity=.0001;
    29         end
     21disp('      creating temperature');
     22md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    3023
    31         %dynamics
    32         md.dt=1*md.yts; %1 year
    33         md.ndt=md.dt*10;
    34         md.artificial_diffusivity=1;
     24disp('      creating flow law paramter');
     25md.B=paterson(md.observed_temperature);
     26md.n=3*ones(md.numberofelements,1);
    3527
    36         %control
    37         md.control_type={'drag'}; %'drag', 'B'
    38         md.nsteps=5;
    39         md.tolx=10^-4;
    40         md.maxiter=20;
    41         md.optscal=10;
    42         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    43         md.meanvel=1000/md.yts; %1000 meters/year
    44         md.epsvel=eps;
     28%Deal with boundary conditions:
     29disp('      boundary conditions for diagnostic model: ');
     30md=SetIceShelfBC(md,'Front.exp');
    4531
    46 
    47         disp('      creating thickness');
    48         hmin=300;
    49         hmax=1000;
    50         ymin=min(md.y);
    51         ymax=max(md.y);
    52         md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
    53         md.firn_layer=10*ones(md.numberofgrids,1);
    54         md.bed=-di*md.thickness+10;
    55         md.surface=md.bed+md.thickness;
    56 
    57         disp('      creating velocities');
    58         md.vx_obs=zeros(md.numberofgrids,1);
    59         md.vy_obs=zeros(md.numberofgrids,1);
    60         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    61        
    62         disp('      creating drag');
    63         md.drag_type=2; %0 none 1 plastic 2 viscous
    64         md.drag=200*ones(md.numberofgrids,1); %q=1.
    65         %Take care of iceshelves: no basal drag
    66         pos=find(md.elementoniceshelf);
    67         md.drag(md.elements(pos,:))=0;
    68         md.p=ones(md.numberofelements,1);
    69         md.q=ones(md.numberofelements,1);
    70 
    71         disp('      creating temperature');
    72         md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    73 
    74         disp('      creating flow law paramter');
    75         md.B=paterson(md.observed_temperature);
    76         md.n=3*ones(md.numberofelements,1);
    77 
    78         disp('      creating accumulation rates');
    79         md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
    80         md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
    81 
    82         %Deal with boundary conditions:
    83        
    84         disp('      boundary conditions for diagnostic model: ');
    85         %Build gridonicefront, array of boundary grids belonging to the icefront:
    86         gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    87         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    88 
    89         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    90         pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
    91         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    92 
    93         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    94         md.segmentonneumann_diag=md.segments(pos,:);
    95         md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    96 
    97         disp('      boundary conditions for prognostic model: ');
    98         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    99         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    100         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    101         md.segmentonneumann_prog=md.segments(pos,:);
    102         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    103         md.neumannvalues_prog(:)=NaN; %free radiation
    104        
    105         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    106         md.segmentonneumann_prog2=md.segments(pos,:);
    107         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    108         md.neumannvalues_prog2(:)=NaN; %free radiation
    109        
    110         disp('      boundary conditions for thermal model: ');
    111         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    112         md.dirichletvalues_thermal=md.observed_temperature;
    113         md.geothermalflux=zeros(md.numberofgrids,1);
    114         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    115 
    116 
    117        
    118 
    119 % Some Cielo code, ignore.
    120 if strcmp(md.cluster,'yes')
    121         ServerDisconnect;
    122 end   
     32%Parallel options
     33md.np=3;
     34md.time=50;
     35md.waitonlock=1;
Note: See TracChangeset for help on using the changeset viewer.