Changeset 904


Ignore:
Timestamp:
06/11/09 14:34:20 (16 years ago)
Author:
seroussi
Message:

simplified Square.par

Location:
issm/trunk/test/Validation/ISMIP
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/test/Validation/ISMIP/TestA/Square.par

    r901 r904  
    2929md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    3030md.segmentonneumann_diag=zeros(0,3);
     31md.gridondirichlet_thermal=zeros(md.numberofgrids,1);

  • issm/trunk/test/Validation/ISMIP/TestB/Square.par

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

  • issm/trunk/test/Validation/ISMIP/TestC/Square.par

    r16 r904  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.8;
    6         md.rho_ice=910;
    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');
     4md.surface=2000-md.x*tan(0.1*pi/180); %to have z>0
     5md.bed=md.surface-1000;
     6md.thickness=md.surface-md.bed;
     7md.firn_layer=0*ones(md.numberofgrids,1);
    138
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     9disp('      creating drag');
     10md.drag_type=2; %0 none 1 plastic 2 viscous
     11md.drag=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)).*sin(md.y*2*pi/max(md.x)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
     12%Take care of iceshelves: no basal drag
     13pos=find(md.elementoniceshelf);
     14md.drag(md.elements(pos,:))=0;
     15md.p=ones(md.numberofelements,1);
     16md.q=ones(md.numberofelements,1);
    2017
    21         %statics
    22         md.lowmem=1;
    23         md.eps_abs=NaN;
    24         md.eps_rel=0.11;
    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
     18disp('      creating flow law paramter');
     19md.B=6.8067*10^7*ones(md.numberofgrids,1);
     20md.n=3*ones(md.numberofelements,1);
    3221
    33         %dynamics
    34         md.dt=1*md.yts; %1 year
    35         md.ndt=md.dt*10;
    36         md.artificial_diffusivity=1;
    37         md.viscosity_overshoot=0;
    38 
    39         %control
    40         md.control_type={'drag'}; %'drag', 'B'
    41         md.nsteps=5;
    42         md.tolx=10^-4;
    43         md.maxiter=20;
    44         md.optscal=10;
    45         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    46         md.meanvel=1000/md.yts; %1000 meters/year
    47         md.epsvel=eps;
    48 
    49 
    50         disp('      creating thickness');
    51         md.surface=2000-md.x*tan(0.1*pi/180); %to have z>0
    52         md.bed=md.surface-1000;
    53         md.thickness=md.surface-md.bed;
    54         md.firn_layer=0*ones(md.numberofgrids,1);
    55 
    56         disp('      creating velocities');
    57         md.vx_obs=zeros(md.numberofgrids,1);
    58         md.vy_obs=zeros(md.numberofgrids,1);
    59         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    60        
    61         disp('      creating drag');
    62         md.drag_type=2; %0 none 1 plastic 2 viscous
    63         md.drag=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)).*sin(md.y*2*pi/max(md.x)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
    64 
    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=6.8067*10^7*ones(md.numberofgrids,1);
    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         %Create grid on boundary fist (because wi can not use mesh)
    85         md.gridonboundary=zeros(md.numberofgrids,1);
    86         pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
    87         md.gridonboundary(pos)=1;
    88 
    89         disp('      boundary conditions for diagnostic model: ');
    90         %Build gridonicefront, array of boundary grids belonging to the icefront:
    91         gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    92         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    93 
    94         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    95         pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
    96         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    97 
    98         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    99         %md.segmentonneumann_diag=md.segments(pos,:);
    100         md.segmentonneumann_diag=zeros(0,3);
    101         %md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    102 
    103         disp('      boundary conditions for prognostic model: ');
    104         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    105         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    106         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    107         %md.segmentonneumann_prog=md.segments(pos,:);
    108         md.segmentonneumann_prog=zeros(0,3);
    109         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    110         md.neumannvalues_prog(:)=NaN; %free radiation
    111        
    112         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    113         %md.segmentonneumann_prog2=md.segments(pos,:);
    114         md.segmentonneumann_prog2=zeros(0,3);
    115         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    116         md.neumannvalues_prog2(:)=NaN; %free radiation
    117        
    118         disp('      boundary conditions for thermal model: ');
    119         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    120         md.dirichletvalues_thermal=md.observed_temperature;
    121         md.geothermalflux=zeros(md.numberofgrids,1);
    122         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    123 
    124 
    125        
    126 
    127 % Some Cielo code, ignore.
    128 if strcmp(md.cluster,'yes')
    129         ServerDisconnect;
    130 end   
     22disp('      boundary conditions for diagnostic model: ');
     23%Create grid on boundary fist (because wi can not use mesh)
     24md.gridonboundary=zeros(md.numberofgrids,1);
     25pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
     26md.gridonboundary(pos)=1;
     27md.gridondirichlet_diag=zeros(md.numberofgrids,1);
     28pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
     29md.dirichletvalues_diag=zeros(md.numberofgrids,2);
     30md.segmentonneumann_diag=zeros(0,3);
     31md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature

  • issm/trunk/test/Validation/ISMIP/TestD/Square.par

    r16 r904  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.8;
    6         md.rho_ice=910;
    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');
     4md.surface=2000-md.x*tan(0.1*pi/180); %to have z>0
     5md.bed=md.surface-1000;
     6md.thickness=md.surface-md.bed;
     7md.firn_layer=0*ones(md.numberofgrids,1);
    138
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     9disp('      creating drag');
     10md.drag_type=2; %0 none 1 plastic 2 viscous
     11md.drag=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
     12%Take care of iceshelves: no basal drag
     13pos=find(md.elementoniceshelf);
     14md.drag(md.elements(pos,:))=0;
     15md.p=ones(md.numberofelements,1);
     16md.q=ones(md.numberofelements,1);
    2017
    21         %statics
    22         md.lowmem=1;
    23         md.eps_abs=NaN;
    24         md.eps_rel=0.01;
    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
     18disp('      creating flow law paramter');
     19md.B=6.8067*10^7*ones(md.numberofgrids,1);
     20md.n=3*ones(md.numberofelements,1);
    3221
    33         %dynamics
    34         md.dt=1*md.yts; %1 year
    35         md.ndt=md.dt*10;
    36         md.artificial_diffusivity=1;
    37         md.viscosity_overshoot=0;
    38 
    39         %control
    40         md.control_type={'drag'}; %'drag', 'B'
    41         md.nsteps=5;
    42         md.tolx=10^-4;
    43         md.maxiter=20;
    44         md.optscal=10;
    45         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    46         md.meanvel=1000/md.yts; %1000 meters/year
    47         md.epsvel=eps;
    48 
    49 
    50         disp('      creating thickness');
    51         md.surface=2000-md.x*tan(0.1*pi/180); %to have z>0
    52         md.bed=md.surface-1000;
    53         md.thickness=md.surface-md.bed;
    54         md.firn_layer=0*ones(md.numberofgrids,1);
    55 
    56         disp('      creating velocities');
    57         md.vx_obs=zeros(md.numberofgrids,1);
    58         md.vy_obs=zeros(md.numberofgrids,1);
    59         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    60        
    61         disp('      creating drag');
    62         md.drag_type=2; %0 none 1 plastic 2 viscous
    63         md.drag=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
    64 
    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=6.8067*10^7*ones(md.numberofgrids,1);
    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         %Create grid on boundary fist (because wi can not use mesh)
    85         md.gridonboundary=zeros(md.numberofgrids,1);
    86         pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
    87         md.gridonboundary(pos)=1;
    88 
    89         disp('      boundary conditions for diagnostic model: ');
    90         %Build gridonicefront, array of boundary grids belonging to the icefront:
    91         gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    92         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    93 
    94         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    95         pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
    96         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    97 
    98         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    99         %md.segmentonneumann_diag=md.segments(pos,:);
    100         md.segmentonneumann_diag=zeros(0,3);
    101         %md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    102 
    103         disp('      boundary conditions for prognostic model: ');
    104         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    105         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    106         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    107         %md.segmentonneumann_prog=md.segments(pos,:);
    108         md.segmentonneumann_prog=zeros(0,3);
    109         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    110         md.neumannvalues_prog(:)=NaN; %free radiation
    111        
    112         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    113         %md.segmentonneumann_prog2=md.segments(pos,:);
    114         md.segmentonneumann_prog2=zeros(0,3);
    115         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    116         md.neumannvalues_prog2(:)=NaN; %free radiation
    117        
    118         disp('      boundary conditions for thermal model: ');
    119         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    120         md.dirichletvalues_thermal=md.observed_temperature;
    121         md.geothermalflux=zeros(md.numberofgrids,1);
    122         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    123 
    124 
    125        
    126 
    127 % Some Cielo code, ignore.
    128 if strcmp(md.cluster,'yes')
    129         ServerDisconnect;
    130 end   
     22disp('      boundary conditions for diagnostic model: ');
     23%Create grid on boundary fist (because wi can not use mesh)
     24md.gridonboundary=zeros(md.numberofgrids,1);
     25pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
     26md.gridonboundary(pos)=1;
     27md.gridondirichlet_diag=zeros(md.numberofgrids,1);
     28pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
     29md.dirichletvalues_diag=zeros(md.numberofgrids,2);
     30md.segmentonneumann_diag=zeros(0,3);
     31md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature

  • issm/trunk/test/Validation/ISMIP/TestE/Square.par

    r16 r904  
    1 
    21%Ok, start defining model parameters here
    32
    4 %material parameters
    5         md.g=9.8;
    6         md.rho_ice=910;
    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');
     4data=load('data.mat');
     5data=data.data;
     6md.surface=zeros(md.numberofgrids,1);
     7md.bed=zeros(md.numberofgrids,1);
     8for i=1:md.numberofgrids
     9        y=md.y(i);
     10        point1=floor(y/100)+1;
     11        point2=min(point1+1,51);
     12        coeff=(y-(point1-1)*100)/100;
     13        md.bed(i)=(1-coeff)*data(point1,2)+coeff*data(point2,2);
     14        md.surface(i)=(1-coeff)*data(point1,3)+coeff*data(point2,3);
     15end
     16md.thickness=md.surface-md.bed;
     17md.firn_layer=0*ones(md.numberofgrids,1);
     18md.thickness(find(~md.thickness))=0.01;
     19md.bed=md.surface-md.thickness;
    1320
    14 %Solution parameters
    15         %parallelization
    16         md.cluster='none';
    17         md.np=2;
    18         md.time=1;
    19         md.exclusive=0;
     21disp('      creating drag');
     22md.drag_type=2; %0 none 1 plastic 2 viscous
     23md.drag=zeros(md.numberofgrids,1);
     24%Take care of iceshelves: no basal drag
     25pos=find(md.elementoniceshelf);
     26md.drag(md.elements(pos,:))=0;
     27md.p=ones(md.numberofelements,1);
     28md.q=ones(md.numberofelements,1);
    2029
    21         %statics
    22         md.lowmem=1;
    23         md.eps_abs=NaN;
    24         md.eps_rel=0.01;
    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
     30disp('      creating flow law paramter');
     31md.B=6.8067*10^7*ones(md.numberofgrids,1);
     32md.n=3*ones(md.numberofelements,1);
    3233
    33         %dynamics
    34         md.dt=1*md.yts; %1 year
    35         md.ndt=md.dt*10;
    36         md.artificial_diffusivity=1;
    37         md.viscosity_overshoot=0;
    38 
    39         %control
    40         md.control_type={'drag'}; %'drag', 'B'
    41         md.nsteps=5;
    42         md.tolx=10^-4;
    43         md.maxiter=20;
    44         md.optscal=10;
    45         md.fit='logarithmic'; %'absolute','relative','logarithmic'
    46         md.meanvel=1000/md.yts; %1000 meters/year
    47         md.epsvel=eps;
    48 
    49 
    50         disp('      creating thickness');
    51         data=load('data.mat');
    52         data=data.data;
    53         md.surface=zeros(md.numberofgrids,1);
    54         md.bed=zeros(md.numberofgrids,1);
    55         for i=1:md.numberofgrids
    56                 y=md.y(i);
    57                 point1=floor(y/100)+1;
    58                 point2=min(point1+1,51);
    59                 coeff=(y-(point1-1)*100)/100;
    60                 md.bed(i)=(1-coeff)*data(point1,2)+coeff*data(point2,2);
    61                 md.surface(i)=(1-coeff)*data(point1,3)+coeff*data(point2,3);
    62         end
    63         md.thickness=md.surface-md.bed;
    64         md.firn_layer=0*ones(md.numberofgrids,1);
    65         md.thickness(find(~md.thickness))=0.01;
    66         md.bed=md.surface-md.thickness;
    67 
    68         disp('      creating velocities');
    69         md.vx_obs=zeros(md.numberofgrids,1);
    70         md.vy_obs=zeros(md.numberofgrids,1);
    71         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    72        
    73         disp('      creating drag');
    74         md.drag_type=2; %0 none 1 plastic 2 viscous
    75         md.drag=zeros(md.numberofgrids,1);
    76         %Take care of iceshelves: no basal drag
    77         pos=find(md.elementoniceshelf);
    78         md.drag(md.elements(pos,:))=0;
    79         md.p=ones(md.numberofelements,1);
    80         md.q=ones(md.numberofelements,1);
    81 
    82         disp('      creating temperature');
    83         md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
    84 
    85         disp('      creating flow law paramter');
    86         md.B=6.8067*10^7*ones(md.numberofgrids,1);
    87         md.n=3*ones(md.numberofelements,1);
    88 
    89         disp('      creating accumulation rates');
    90         md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
    91         md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
    92 
    93         %Deal with boundary conditions:
    94        
    95         %Create grid on boundary fist (because wi can not use mesh)
    96         md.gridonboundary=zeros(md.numberofgrids,1);
    97         pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
    98         md.gridonboundary(pos)=1;
    99 
    100         disp('      boundary conditions for diagnostic model: ');
    101         %Build gridonicefront, array of boundary grids belonging to the icefront:
    102         gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
    103         gridonicefront=double(md.gridonboundary & gridinsideicefront);
    104 
    105         md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    106         pos=find(md.y==0 | md.y==max(md.y));md.gridondirichlet_diag(pos)=1;
    107         md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    108 
    109         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    110         %md.segmentonneumann_diag=md.segments(pos,:);
    111         md.segmentonneumann_diag=zeros(0,3);
    112         %md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
    113 
    114         disp('      boundary conditions for prognostic model: ');
    115         md.gridondirichlet_prog=zeros(md.numberofgrids,1);
    116         md.dirichletvalues_prog=zeros(md.numberofgrids,1);
    117         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    118         %md.segmentonneumann_prog=md.segments(pos,:);
    119         md.segmentonneumann_prog=zeros(0,3);
    120         md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
    121         md.neumannvalues_prog(:)=NaN; %free radiation
    122        
    123         %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
    124         %md.segmentonneumann_prog2=md.segments(pos,:);
    125         md.segmentonneumann_prog2=zeros(0,3);
    126         md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
    127         md.neumannvalues_prog2(:)=NaN; %free radiation
    128        
    129         disp('      boundary conditions for thermal model: ');
    130         md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    131         md.dirichletvalues_thermal=md.observed_temperature;
    132         md.geothermalflux=zeros(md.numberofgrids,1);
    133         pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
    134 
    135 
    136        
    137 
    138 % Some Cielo code, ignore.
    139 if strcmp(md.cluster,'yes')
    140         ServerDisconnect;
    141 end   
     34disp('      boundary conditions for diagnostic model: ');
     35%Create grid on boundary fist (because wi can not use mesh)
     36md.gridonboundary=zeros(md.numberofgrids,1);
     37pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
     38md.gridonboundary(pos)=1;
     39md.gridondirichlet_diag=zeros(md.numberofgrids,1);
     40pos=find(md.y==0 | md.y==max(md.y));md.gridondirichlet_diag(pos)=1;
     41md.dirichletvalues_diag=zeros(md.numberofgrids,2);
     42md.segmentonneumann_diag=zeros(0,3);
     43md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
Note: See TracChangeset for help on using the changeset viewer.