Changeset 900


Ignore:
Timestamp:
06/11/09 13:19:18 (15 years ago)
Author:
Mathieu Morlighem
Message:

simplified Square.par

Location:
issm/trunk/test/Validation
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/test/Validation/EISMINT/Diagnostic1/Square.par

    r899 r900  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.81;
    6         md.rho_ice=917;
    7         md.rho_water=1028;
    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');
     4ymin=min(md.y);
     5ymax=max(md.y);
     6md.thickness=500*ones(md.numberofgrids,1);
     7md.firn_layer=0*ones(md.numberofgrids,1);
     8md.bed=-di*md.thickness;
     9md.surface=md.bed+md.thickness;
    1310
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     11disp('      creating drag');
     12md.drag_type=2; %0 none 1 plastic 2 viscous
     13md.drag=200*ones(md.numberofgrids,1); %q=1.
     14%Take care of iceshelves: no basal drag
     15pos=find(md.elementoniceshelf);
     16md.drag(md.elements(pos,:))=0;
     17md.p=ones(md.numberofelements,1);
     18md.q=ones(md.numberofelements,1);
    2019
    21         %statics
    22         md.lowmem=1;
    23         md.eps_rel=0.01;
    24         md.eps_abs=10;
    25         md.penalty_offset=3;
    26         md.penalty_melting=10^7;
    27         if md.numberofgrids<1000000,
    28         md.sparsity=.001;
    29         else
    30         md.sparsity=.0001;
    31         end
     20disp('      creating temperature');
     21md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    3222
    33         %dynamics
    34         md.dt=1*md.yts; %1 year
    35         md.ndt=md.dt*10;
    36         md.artificial_diffusivity=1;
     23disp('      creating flow law paramter');
     24md.B=1.7687*10^8*ones(md.numberofgrids,1);
     25md.n=3*ones(md.numberofelements,1);
    3726
    38         %control
    39         md.control_type={'drag'}; %'drag', 'B'
    40         md.nsteps=5;
    41         md.tolx=10^-4;
    42         md.maxiter=20;
    43         md.optscal=10;
    44         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    45         md.meanvel=1000/md.yts; %1000 meters/year
    46         md.epsvel=eps;
     27disp('      creating accumulation rates');
     28md.accumulation=0.2*ones(md.numberofgrids,1); %0.2m/a
     29md.melting=0*ones(md.numberofgrids,1); %0m/a
    4730
    48 
    49         disp('      creating thickness');
    50         ymin=min(md.y);
    51         ymax=max(md.y);
    52         md.thickness=500*ones(md.numberofgrids,1);
    53         md.firn_layer=0*ones(md.numberofgrids,1);
    54         md.bed=-di*md.thickness;
    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.B=1.7687*10^8*ones(md.numberofgrids,1);
    77         md.n=3*ones(md.numberofelements,1);
    78 
    79         disp('      creating accumulation rates');
    80         md.accumulation=0.2*ones(md.numberofgrids,1); %10m/a
    81         md.melting=0*ones(md.numberofgrids,1); %1m/a
    82 
    83         %Deal with boundary conditions:
    84        
    85         disp('      boundary conditions for diagnostic model: ');
    86         %Build gridonicefront, array of boundary grids belonging to the icefront:
    87         gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    88         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    89 
    90         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    91         pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
    92         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    93 
    94         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    95         md.segmentonneumann_diag=md.segments(pos,:);
    96         md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    97 
    98         disp('      boundary conditions for prognostic model: ');
    99         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    100         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    101         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    102         md.segmentonneumann_prog=md.segments(pos,:);
    103         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    104         md.neumannvalues_prog(:)=NaN; %free radiation
    105        
    106         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    107         md.segmentonneumann_prog2=md.segments(pos,:);
    108         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    109         md.neumannvalues_prog2(:)=NaN; %free radiation
    110        
    111         disp('      boundary conditions for thermal model: ');
    112         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    113         md.dirichletvalues_thermal=md.observed_temperature;
    114         md.geothermalflux=zeros(md.numberofgrids,1);
    115         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    116 
    117 
    118        
    119 
    120 % Some Cielo code, ignore.
    121 if strcmp(md.cluster,'yes')
    122         ServerDisconnect;
    123 end   
     31disp('      boundary conditions ');
     32md=SetMarineIceSheetBC(md,'Front.exp');

  • issm/trunk/test/Validation/EISMINT/Diagnostic2/Square.par

    r899 r900  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.81;
    6         md.rho_ice=917;
    7         md.rho_water=1028;
    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');
     4ymin=min(md.y);
     5ymax=max(md.y);
     6md.thickness=500*ones(md.numberofgrids,1);
     7md.firn_layer=0*ones(md.numberofgrids,1);
     8md.bed=-di*md.thickness;
     9md.surface=md.bed+md.thickness;
    1310
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     11disp('      creating drag');
     12md.drag_type=2; %0 none 1 plastic 2 viscous
     13md.drag=200*ones(md.numberofgrids,1); %q=1.
     14%Take care of iceshelves: no basal drag
     15pos=find(md.elementoniceshelf);
     16md.drag(md.elements(pos,:))=0;
     17md.p=ones(md.numberofelements,1);
     18md.q=ones(md.numberofelements,1);
    2019
    21         %statics
    22         md.lowmem=1;
    23         md.eps_rel=0.01;
    24         md.eps_abs=10;
    25         md.penalty_offset=3;
    26         md.penalty_melting=10^7;
    27         if md.numberofgrids<1000000,
    28         md.sparsity=.001;
    29         else
    30         md.sparsity=.0001;
    31         end
     20disp('      creating temperature');
     21md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    3222
    33         %dynamics
    34         md.dt=1*md.yts; %1 year
    35         md.ndt=md.dt*10;
    36         md.artificial_diffusivity=1;
     23disp('      creating flow law paramter');
     24md.B=1.7687*10^8*ones(md.numberofgrids,1);
     25md.n=3*ones(md.numberofelements,1);
    3726
    38         %control
    39         md.control_type={'drag'}; %'drag', 'B'
    40         md.nsteps=5;
    41         md.tolx=10^-4;
    42         md.maxiter=20;
    43         md.optscal=10;
    44         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    45         md.meanvel=1000/md.yts; %1000 meters/year
    46         md.epsvel=eps;
     27disp('      creating accumulation rates');
     28md.accumulation=0.2*ones(md.numberofgrids,1); %0.2m/a
     29md.melting=0*ones(md.numberofgrids,1); %0m/a
    4730
    48 
    49         disp('      creating thickness');
    50         ymin=min(md.y);
    51         ymax=max(md.y);
    52         md.thickness=500*ones(md.numberofgrids,1);
    53         md.firn_layer=0*ones(md.numberofgrids,1);
    54         md.bed=-di*md.thickness;
    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.B=1.7687*10^8*ones(md.numberofgrids,1);
    77         md.n=3*ones(md.numberofelements,1);
    78 
    79         disp('      creating accumulation rates');
    80         md.accumulation=0.2*ones(md.numberofgrids,1); %10m/a
    81         md.melting=0*ones(md.numberofgrids,1); %1m/a
    82 
    83         %Deal with boundary conditions:
    84        
    85         disp('      boundary conditions for diagnostic model: ');
    86         %Build gridonicefront, array of boundary grids belonging to the icefront:
    87         gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    88         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    89 
    90         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    91         pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
    92         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    93 
    94         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    95         md.segmentonneumann_diag=md.segments(pos,:);
    96         md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    97 
    98         disp('      boundary conditions for prognostic model: ');
    99         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    100         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    101         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    102         pos=find(md.y==max(md.y)); %grids on the upper boundary condition
    103         md.dirichletvalues_diag(pos,2)=400*(((md.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);
    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   
     31disp('      boundary conditions ');
     32md=SetMarineIceSheetBC(md,'Front.exp');
     33pos=find(md.y==max(md.y)); %grids on the upper boundary condition
     34md.dirichletvalues_diag(pos,2)=400*(((md.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);

  • issm/trunk/test/Validation/EISMINT/MassConservation/Square.par

    r899 r900  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.81;
    6         md.rho_ice=917;
    7         md.rho_water=1028;
    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');
     4ymin=min(md.y);
     5ymax=max(md.y);
     6md.thickness=500*ones(md.numberofgrids,1);
     7md.firn_layer=0*ones(md.numberofgrids,1);
     8md.bed=-di*md.thickness;
     9md.surface=md.bed+md.thickness;
    1310
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     11disp('      creating drag');
     12md.drag_type=2; %0 none 1 plastic 2 viscous
     13md.drag=200*ones(md.numberofgrids,1); %q=1.
     14%Take care of iceshelves: no basal drag
     15pos=find(md.elementoniceshelf);
     16md.drag(md.elements(pos,:))=0;
     17md.p=ones(md.numberofelements,1);
     18md.q=ones(md.numberofelements,1);
    2019
    21         %statics
    22         md.lowmem=1;
    23         md.eps_rel=0.01;
    24         md.eps_abs=10;
    25         md.penalty_offset=3;
    26         md.penalty_melting=10^7;
    27         if md.numberofgrids<1000000,
    28         md.sparsity=.001;
    29         else
    30         md.sparsity=.0001;
    31         end
     20disp('      creating temperature');
     21md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    3222
    33         %dynamics
    34         md.dt=1*md.yts; %1 year
    35         md.ndt=md.dt*10;
    36         md.artificial_diffusivity=1;
     23disp('      creating flow law paramter');
     24%md.B=paterson(md.observed_temperature);
     25md.B=1.7687*10^8*ones(md.numberofgrids,1);
     26md.n=3*ones(md.numberofelements,1);
    3727
    38         %control
    39         md.control_type={'drag'}; %'drag', 'B'
    40         md.nsteps=5;
    41         md.tolx=10^-4;
    42         md.maxiter=20;
    43         md.optscal=10;
    44         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    45         md.meanvel=1000/md.yts; %1000 meters/year
    46         md.epsvel=eps;
     28disp('      creating accumulation rates');
     29md.accumulation=0.0*ones(md.numberofgrids,1); %0m/a
     30md.melting=0*ones(md.numberofgrids,1); %0m/a
    4731
     32disp('      boundary conditions ');
     33md=SetMarineIceSheetBC(md,'Front.exp');
    4834
    49         disp('      creating thickness');
    50         ymin=min(md.y);
    51         ymax=max(md.y);
    52         md.thickness=500*ones(md.numberofgrids,1);
    53         md.firn_layer=0*ones(md.numberofgrids,1);
    54         md.bed=-di*md.thickness;
    55         md.surface=md.bed+md.thickness;
     35%analytical test
     36pos=find(md.y==200000); %grids on the upper boundary condition
     37md.gridondirichlet_diag=ones(md.numberofgrids,1);
     38md.dirichletvalues_diag(:,1)=0;
     39md.dirichletvalues_diag(:,2)=-400;
    5640
    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.B=1.7687*10^8*ones(md.numberofgrids,1);
    77         md.n=3*ones(md.numberofelements,1);
    78 
    79         disp('      creating accumulation rates');
    80         md.accumulation=0.0*ones(md.numberofgrids,1); %10m/a
    81         md.melting=0*ones(md.numberofgrids,1); %1m/a
    82 
    83         %Deal with boundary conditions:
    84        
    85         disp('      boundary conditions for diagnostic model: ');
    86         %Build gridonicefront, array of boundary grids belonging to the icefront:
    87         gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    88         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    89 
    90         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    91         pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
    92         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    93 
    94         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    95         md.segmentonneumann_diag=md.segments(pos,:);
    96         md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    97 
    98         %analytical test
    99         pos=find(md.y==200000); %grids on the upper boundary condition
    100         md.gridondirichlet_diag=ones(md.numberofgrids,1);
    101         md.dirichletvalues_diag(:,1)=0;
    102         md.dirichletvalues_diag(:,2)=-400;
    103 
    104         disp('      boundary conditions for prognostic model: ');
    105         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    106         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    107         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    108         md.segmentonneumann_prog=md.segments(pos,:);
    109         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    110         md.neumannvalues_prog(:)=NaN; %free radiation
    111        
    112         %analytical test
    113         pos=find(md.y==200000); %grids on the upper boundary condition
    114         md.gridondirichlet_prog(pos)=1;
    115         md.dirichletvalues_prog(pos)=500;
    116 
    117         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    118         md.segmentonneumann_prog2=md.segments(pos,:);
    119         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    120         md.neumannvalues_prog2(:)=NaN; %free radiation
    121        
    122         disp('      boundary conditions for thermal model: ');
    123         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    124         md.dirichletvalues_thermal=md.observed_temperature;
    125         md.geothermalflux=zeros(md.numberofgrids,1);
    126         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    127 
    128 
    129        
    130 
    131 % Some Cielo code, ignore.
    132 if strcmp(md.cluster,'yes')
    133         ServerDisconnect;
    134 end   
     41%analytical test
     42pos=find(md.y==200000); %grids on the upper boundary condition
     43md.gridondirichlet_prog(pos)=1;
     44md.dirichletvalues_prog(pos)=500;

  • issm/trunk/test/Validation/EISMINT/Transient/Square.par

    r899 r900  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.81;
    6         md.rho_ice=917;
    7         md.rho_water=1028;
    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');
     4ymin=min(md.y);
     5ymax=max(md.y);
     6md.thickness=500*ones(md.numberofgrids,1);
     7md.firn_layer=0*ones(md.numberofgrids,1);
     8md.bed=-di*md.thickness;
     9md.surface=md.bed+md.thickness;
    1310
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     11disp('      creating drag');
     12md.drag_type=2; %0 none 1 plastic 2 viscous
     13md.drag=200*ones(md.numberofgrids,1); %q=1.
     14%Take care of iceshelves: no basal drag
     15pos=find(md.elementoniceshelf);
     16md.drag(md.elements(pos,:))=0;
     17md.p=ones(md.numberofelements,1);
     18md.q=ones(md.numberofelements,1);
    2019
    21         %statics
    22         md.lowmem=1;
    23         md.eps_rel=0.01;
    24         md.eps_abs=10;
    25         md.penalty_offset=3;
    26         md.penalty_melting=10^7;
    27         if md.numberofgrids<1000000,
    28         md.sparsity=.001;
    29         else
    30         md.sparsity=.0001;
    31         end
     20disp('      creating temperature');
     21md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    3222
    33         %dynamics
    34         md.dt=1*md.yts; %1 year
    35         md.ndt=md.dt*10;
    36         md.artificial_diffusivity=1;
     23disp('      creating flow law paramter');
     24md.B=1.7687*10^8*ones(md.numberofgrids,1);
     25md.n=3*ones(md.numberofelements,1);
    3726
    38         %control
    39         md.control_type={'drag'}; %'drag', 'B'
    40         md.nsteps=5;
    41         md.tolx=10^-4;
    42         md.maxiter=20;
    43         md.optscal=10;
    44         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    45         md.meanvel=1000/md.yts; %1000 meters/year
    46         md.epsvel=eps;
     27disp('      creating accumulation rates');
     28md.accumulation=0.2*ones(md.numberofgrids,1); %0.2m/a
     29md.melting=0*ones(md.numberofgrids,1); %0m/a
    4730
     31disp('      boundary conditions ');
     32md=SetMarineIceSheetBC(md,'Front.exp');
    4833
    49         disp('      creating thickness');
    50         ymin=min(md.y);
    51         ymax=max(md.y);
    52         md.thickness=500*ones(md.numberofgrids,1);
    53         md.firn_layer=0*ones(md.numberofgrids,1);
    54         md.bed=-di*md.thickness;
    55         md.surface=md.bed+md.thickness;
     34disp('      boundary conditions for prognostic model: ');
     35md.gridondirichlet_prog=zeros(md.numberofgrids,1);
     36md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    5637
    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);
     38pos=find(md.y==max(md.y)); %grids on the upper boundary condition
     39md.dirichletvalues_diag(pos,2)=400*(((md.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);
     40md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
     41md.neumannvalues_prog(:)=NaN; %free radiation
    7042
    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.B=1.7687*10^8*ones(md.numberofgrids,1);
    77         md.n=3*ones(md.numberofelements,1);
    78 
    79         disp('      creating accumulation rates');
    80         md.accumulation=0.2*ones(md.numberofgrids,1); %10m/a
    81         md.melting=0*ones(md.numberofgrids,1); %1m/a
    82 
    83         %Deal with boundary conditions:
    84        
    85         disp('      boundary conditions for diagnostic model: ');
    86         %Build gridonicefront, array of boundary grids belonging to the icefront:
    87         gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    88         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    89 
    90         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    91         pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
    92         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    93 
    94         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    95         md.segmentonneumann_diag=md.segments(pos,:);
    96         md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    97 
    98         disp('      boundary conditions for prognostic model: ');
    99         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    100         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    101         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    102         pos=find(md.y==max(md.y)); %grids on the upper boundary condition
    103         md.dirichletvalues_diag(pos,2)=400*(((md.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);
    104         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    105         md.neumannvalues_prog(:)=NaN; %free radiation
    106        
    107         pos=find(md.y==max(md.y)); %grids on the upper boundary condition
    108         md.gridondirichlet_prog(pos)=1;
    109         md.dirichletvalues_prog(pos)=500;
    110 
    111         pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    112         md.segmentonneumann_prog2=md.segments(pos,:);
    113         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    114         md.neumannvalues_prog2(:)=NaN; %free radiation
    115        
    116         disp('      boundary conditions for thermal model: ');
    117         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    118         md.dirichletvalues_thermal=md.observed_temperature;
    119         md.geothermalflux=zeros(md.numberofgrids,1);
    120         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    121 
    122 
    123        
    124 
    125 % Some Cielo code, ignore.
    126 if strcmp(md.cluster,'yes')
    127         ServerDisconnect;
    128 end   
     43pos=find(md.y==max(md.y)); %grids on the upper boundary condition
     44md.gridondirichlet_prog(pos)=1;
     45md.dirichletvalues_prog(pos)=500;

  • issm/trunk/test/Validation/ThermalTests/Melting/Square.par

    r544 r900  
    11%Ok, start defining model parameters here
    22
    3         %parallelization
    4         md.cluster='none';
     3%parallelization
     4md.cluster='none';
    55
    6         disp('      creating thickness');
    7         h=1000;
    8         md.thickness=h*ones(md.numberofgrids,1);
    9         md.firn_layer=10*ones(md.numberofgrids,1);
    10         md.bed=-1000*ones(md.numberofgrids,1);
    11         md.surface=md.bed+md.thickness;
     6disp('      creating thickness');
     7h=1000;
     8md.thickness=h*ones(md.numberofgrids,1);
     9md.firn_layer=10*ones(md.numberofgrids,1);
     10md.bed=-1000*ones(md.numberofgrids,1);
     11md.surface=md.bed+md.thickness;
    1212
    13         disp('      creating velocities');
    14         md.vx_obs=zeros(md.numberofgrids,1);
    15         md.vy_obs=zeros(md.numberofgrids,1);
    16         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    17        
    18         disp('      creating drag');
    19         md.drag_type=2; %0 none 1 plastic 2 viscous
    20         md.drag=200*ones(md.numberofgrids,1); %q=1.
    21         %Take care of iceshelves: no basal drag
    22         pos=find(md.elementoniceshelf);
    23         md.drag(md.elements(pos,:))=0;
    24         md.p=ones(md.numberofelements,1);
    25         md.q=ones(md.numberofelements,1);
     13disp('      creating velocities');
     14md.vx_obs=zeros(md.numberofgrids,1);
     15md.vy_obs=zeros(md.numberofgrids,1);
     16md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    2617
    27         disp('      creating temperatures');
    28         md.observed_temperature=273.15*ones(md.numberofgrids,1);
     18disp('      creating drag');
     19md.drag_type=2; %0 none 1 plastic 2 viscous
     20md.drag=200*ones(md.numberofgrids,1); %q=1.
     21%Take care of iceshelves: no basal drag
     22pos=find(md.elementoniceshelf);
     23md.drag(md.elements(pos,:))=0;
     24md.p=ones(md.numberofelements,1);
     25md.q=ones(md.numberofelements,1);
    2926
    30         disp('      creating flow law paramter');
    31         md.B=paterson(md.observed_temperature);
    32         md.n=3*ones(md.numberofelements,1);
    33        
    34         disp('      creating accumulation rates');
    35         md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
    36         md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
     27disp('      creating temperatures');
     28md.observed_temperature=273.15*ones(md.numberofgrids,1);
    3729
    38         %Deal with boundary conditions:
    39        
    40         disp('      boundary conditions for diagnostic model');
    41         md=SetIceShelfBC(md,'Front.exp');
    42        
    43         disp('      boundary conditions for thermal model');
    44         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    45         md.dirichletvalues_thermal=md.observed_temperature;
    46         md.geothermalflux=zeros(md.numberofgrids,1);
    47         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=1*10^-3; %1 mW/m^2
     30disp('      creating flow law paramter');
     31md.B=paterson(md.observed_temperature);
     32md.n=3*ones(md.numberofelements,1);
     33
     34disp('      creating accumulation rates');
     35md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
     36md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
     37
     38%Deal with boundary conditions:
     39
     40disp('      boundary conditions for diagnostic model');
     41md=SetIceShelfBC(md,'Front.exp');
     42
     43disp('      boundary conditions for thermal model');
     44md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
     45md.dirichletvalues_thermal=md.observed_temperature;
     46md.geothermalflux=zeros(md.numberofgrids,1);
     47pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=1*10^-3; %1 mW/m^2

  • issm/trunk/test/Validation/ThermalTests/Simpleadvection/Square.par

    r514 r900  
    11%Ok, start defining model parameters here
    22
    3         %parallelization
    4         md.cluster='none';
     3%parallelization
     4md.cluster='none';
    55
    6         disp('      creating thickness');
    7         h=1000;
    8         md.thickness=h*ones(md.numberofgrids,1);
    9         md.firn_layer=10*ones(md.numberofgrids,1);
    10         md.bed=-1000*ones(md.numberofgrids,1);
    11         md.surface=md.bed+md.thickness;
     6disp('      creating thickness');
     7h=1000;
     8md.thickness=h*ones(md.numberofgrids,1);
     9md.firn_layer=10*ones(md.numberofgrids,1);
     10md.bed=-1000*ones(md.numberofgrids,1);
     11md.surface=md.bed+md.thickness;
    1212
    13         disp('      creating velocities');
    14         md.vx_obs=zeros(md.numberofgrids,1);
    15         md.vy_obs=zeros(md.numberofgrids,1);
    16         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    17        
    18         disp('      creating drag');
    19         md.drag_type=2; %0 none 1 plastic 2 viscous
    20         md.drag=200*ones(md.numberofgrids,1); %q=1.
    21         %Take care of iceshelves: no basal drag
    22         pos=find(md.elementoniceshelf);
    23         md.drag(md.elements(pos,:))=0;
    24         md.p=ones(md.numberofelements,1);
    25         md.q=ones(md.numberofelements,1);
     13disp('      creating velocities');
     14md.vx_obs=zeros(md.numberofgrids,1);
     15md.vy_obs=zeros(md.numberofgrids,1);
     16md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    2617
    27         disp('      creating temperatures');
    28         md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
     18disp('      creating drag');
     19md.drag_type=2; %0 none 1 plastic 2 viscous
     20md.drag=200*ones(md.numberofgrids,1); %q=1.
     21%Take care of iceshelves: no basal drag
     22pos=find(md.elementoniceshelf);
     23md.drag(md.elements(pos,:))=0;
     24md.p=ones(md.numberofelements,1);
     25md.q=ones(md.numberofelements,1);
    2926
    30         disp('      creating flow law paramter');
    31         md.B=paterson(md.observed_temperature);
    32         md.n=3*ones(md.numberofelements,1);
    33        
    34         disp('      creating accumulation rates');
    35         md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
    36         md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
     27disp('      creating temperatures');
     28md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    3729
    38         %Deal with boundary conditions:
    39        
    40         disp('      boundary conditions for diagnostic model');
    41         md=SetIceShelfBC(md,'Front.exp');
    42        
    43         disp('      boundary conditions for thermal model');
    44         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    45         md.dirichletvalues_thermal=md.observed_temperature;
    46         md.geothermalflux=zeros(md.numberofgrids,1);
    47         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
     30disp('      creating flow law paramter');
     31md.B=paterson(md.observed_temperature);
     32md.n=3*ones(md.numberofelements,1);
     33
     34disp('      creating accumulation rates');
     35md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
     36md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
     37
     38%Deal with boundary conditions:
     39
     40disp('      boundary conditions for diagnostic model');
     41md=SetIceShelfBC(md,'Front.exp');
     42
     43disp('      boundary conditions for thermal model');
     44md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
     45md.dirichletvalues_thermal=md.observed_temperature;
     46md.geothermalflux=zeros(md.numberofgrids,1);
     47pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2

  • issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux/Square.par

    r515 r900  
    11%Ok, start defining model parameters here
    22
    3         %parallelization
    4         md.cluster='none';
     3%parallelization
     4md.cluster='none';
    55
    6         disp('      creating thickness');
    7         h=1000;
    8         md.thickness=h*ones(md.numberofgrids,1);
    9         md.firn_layer=10*ones(md.numberofgrids,1);
    10         md.bed=zeros(md.numberofgrids,1);
    11         md.surface=md.thickness;
     6disp('      creating thickness');
     7h=1000;
     8md.thickness=h*ones(md.numberofgrids,1);
     9md.firn_layer=10*ones(md.numberofgrids,1);
     10md.bed=zeros(md.numberofgrids,1);
     11md.surface=md.thickness;
    1212
    13         disp('      creating velocities');
    14         md.vx_obs=zeros(md.numberofgrids,1);
    15         md.vy_obs=zeros(md.numberofgrids,1);
    16         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    17        
    18         disp('      creating drag');
    19         md.drag_type=2; %0 none 1 plastic 2 viscous
    20         md.drag=200*ones(md.numberofgrids,1); %q=1.
    21         %Take care of iceshelves: no basal drag
    22         pos=find(md.elementoniceshelf);
    23         md.drag(md.elements(pos,:))=0;
    24         md.p=ones(md.numberofelements,1);
    25         md.q=ones(md.numberofelements,1);
     13disp('      creating velocities');
     14md.vx_obs=zeros(md.numberofgrids,1);
     15md.vy_obs=zeros(md.numberofgrids,1);
     16md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    2617
    27         disp('      creating temperatures');
    28         md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
     18disp('      creating drag');
     19md.drag_type=2; %0 none 1 plastic 2 viscous
     20md.drag=200*ones(md.numberofgrids,1); %q=1.
     21%Take care of iceshelves: no basal drag
     22pos=find(md.elementoniceshelf);
     23md.drag(md.elements(pos,:))=0;
     24md.p=ones(md.numberofelements,1);
     25md.q=ones(md.numberofelements,1);
    2926
    30         disp('      creating flow law paramter');
    31         md.B=paterson(md.observed_temperature);
    32         md.n=3*ones(md.numberofelements,1);
    33        
    34         disp('      creating accumulation rates');
    35         md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
    36         md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
     27disp('      creating temperatures');
     28md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    3729
    38         %Deal with boundary conditions:
    39        
    40         disp('      boundary conditions for diagnostic model');
    41         md=SetIceShelfBC(md,'Front.exp');
    42        
    43         disp('      boundary conditions for thermal model');
    44         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    45         md.dirichletvalues_thermal=md.observed_temperature;
    46         md.geothermalflux=zeros(md.numberofgrids,1);
    47         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=0.1; %100 mW/m^2
     30disp('      creating flow law paramter');
     31md.B=paterson(md.observed_temperature);
     32md.n=3*ones(md.numberofelements,1);
     33
     34disp('      creating accumulation rates');
     35md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
     36md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
     37
     38%Deal with boundary conditions:
     39
     40disp('      boundary conditions for diagnostic model');
     41md=SetIceShelfBC(md,'Front.exp');
     42
     43disp('      boundary conditions for thermal model');
     44md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
     45md.dirichletvalues_thermal=md.observed_temperature;
     46md.geothermalflux=zeros(md.numberofgrids,1);
     47pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=0.1; %100 mW/m^2
Note: See TracChangeset for help on using the changeset viewer.