Changeset 4876


Ignore:
Timestamp:
07/29/10 10:43:53 (15 years ago)
Author:
seroussi
Message:

Validation tests Thermaltests

Location:
issm/trunk/test/Validation/ThermalTests
Files:
12 edited

Legend:

Unmodified
Added
Removed

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

    r3768 r4876  
    33%parallelization
    44md.cluster='none';
     5md.np=3;
     6md.dt=0;
    57
    68disp('      creating thickness');
    … …  
    3941
    4042disp('      boundary conditions for diagnostic model');
    41 md=SetIceShelfBC(md,'Front.exp');
     43md=SetMarineIceSheetBC(md,'Front.exp');
    4244
    4345disp('      boundary conditions for thermal model');
    44 md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    45 md.dirichletvalues_thermal=md.observed_temperature;
     46md.temperature=md.observed_temperature;
     47md.spctemperature(:,1)=ones(md.numberofgrids,1); %surface temperature
     48md.spctemperature(:,2)=md.observed_temperature;
    4649md.geothermalflux=zeros(md.numberofgrids,1);
    4750pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=1*10^-3; %1 mW/m^2

  • issm/trunk/test/Validation/ThermalTests/Melting/runme.m

    r2553 r4876  
    2222
    2323%modeled  results
    24 %md=solve(md,'analysis_type','thermal','sub_analysis_type','steady');
    25 md=solve(md,'analysis_type','thermal','sub_analysis_type','steady');
     24md=solve(md,'analysis_type',ThermalSolutionEnum);
    2625
    2726%plot results
    28 relative=abs((md.results.thermal.melting-melting)./melting)*100;
    29 relative(find(md.results.thermal.melting==melting))=0;
    30 plotmodel(md,'data','melting','title','modeled melting','data',melting,'title','analytical melting','data',md.results.thermal.melting-melting,'title','absolute error','data',relative,'title','relative error [%]','layer#all',1)
     27comp_melting=zeros(md.numberofgrids,1);
     28comp_melting(md.results.ThermalSolution.MeltingRate.index)=md.results.ThermalSolution.MeltingRate.value;
     29relative=abs((comp_melting-melting)./melting)*100;
     30relative(find(comp_melting==melting))=0;
     31plotmodel(md,'data',comp_melting,'title','modeled melting','data',melting,'title','analytical melting','data',comp_melting-melting,'title','absolute error','data',relative,'title','relative error [%]','layer#all',1)

  • issm/trunk/test/Validation/ThermalTests/Melting/thermalboundarycondition.m

    r16 r4876  
    22%Ok, start defining model BC
    33disp('      boundary conditions for thermal model ');
    4 if strcmp(md.type,'3d')
     4if md.dim==3,
    55%       pos2=find(md.elementonsurface); md.gridondirichlet_thermal(md.elements(pos2,4:6))=1; md.dirichletvalues_thermal(md.elements(pos2,4:6))=0;
    66        md.vx=zeros(md.numberofgrids,1);

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

    r3768 r4876  
    33%parallelization
    44md.cluster='none';
     5md.dt=0;
    56
    67disp('      creating thickness');
    … …  
    3738
    3839%Deal with boundary conditions:
    39 
    4040disp('      boundary conditions for diagnostic model');
    41 md=SetIceShelfBC(md,'Front.exp');
     41md=SetMarineIceSheetBC(md,'Front.exp');
    4242
    4343disp('      boundary conditions for thermal model');
    44 md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    45 md.dirichletvalues_thermal=md.observed_temperature;
     44md.temperature=md.observed_temperature;
     45md.spctemperature(:,1)=ones(md.numberofgrids,1); %surface temperature
     46md.spctemperature(:,2)=md.observed_temperature;
    4647md.geothermalflux=zeros(md.numberofgrids,1);
    4748pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
     49
     50disp('      cluster properties');
     51md.np=3;

  • issm/trunk/test/Validation/ThermalTests/Simpleadvection/runme.m

    r899 r4876  
    2323
    2424%modeled  results
    25 md=solve(md,'analysis_type','thermal','sub_analysis_type','steady','package','ice');
     25md=solve(md,'analysis_type',ThermalSolutionEnum);
    2626
    2727%plot results
    28 relative=abs((md.results.thermal.temperature-md.observed_temperature)./md.observed_temperature)*100;
    29 relative(find(md.results.thermal.temperature==md.observed_temperature))=0;
    30 plotmodel(md,'data',md.results.thermal.temperature,'title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',md.results.thermal.temperature-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)
     28comp_temp=zeros(md.numberofgrids,1);
     29comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value;
     30relative=abs((comp_temp-md.observed_temperature)./md.observed_temperature)*100;
     31relative(find(comp_temp==md.observed_temperature))=0;
     32plotmodel(md,'data',comp_temp,'title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',comp_temp-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)

  • issm/trunk/test/Validation/ThermalTests/Simpleadvection/thermalboundarycondition.m

    r16 r4876  
    22%Ok, start defining model BC
    33disp('      boundary conditions for thermal model ');
    4 if strcmp(md.type,'3d')
    5         pos1=find(md.elementonbed);     md.gridondirichlet_thermal(md.elements(pos1,1:3))=1; md.dirichletvalues_thermal(md.elements(pos1,1:3))=10;
    6         pos2=find(md.elementonsurface); md.gridondirichlet_thermal(md.elements(pos2,4:6))=1; md.dirichletvalues_thermal(md.elements(pos2,4:6))=0;
     4if md.dim==3,
     5        pos1=find(md.elementonbed);     md.spctemperature(md.elements(pos1,1:3),1)=1; md.spctemperature(md.elements(pos1,1:3),2)=10;
     6        pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6),1)=1; md.spctemperature(md.elements(pos2,4:6),2)=0;
    77        md.vx=zeros(md.numberofgrids,1);
    88        md.vy=zeros(md.numberofgrids,1);

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

    r3768 r4876  
    11%Ok, start defining model parameters here
    22
    3         %parallelization
    4         md.cluster='none';
     3md.dt=0;
    54
    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;
     5%parallelization
     6md.np=3;
    127
    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_coefficient=200*ones(md.numberofgrids,1); %q=1.
    21         %Take care of iceshelves: no basal drag
    22         pos=find(md.elementoniceshelf);
    23         md.drag_coefficient(md.elements(pos,:))=0;
    24         md.drag_p=ones(md.numberofelements,1);
    25         md.drag_q=ones(md.numberofelements,1);
     8disp('      creating thickness');
     9h=1000;
     10md.thickness=h*ones(md.numberofgrids,1);
     11md.firn_layer=10*ones(md.numberofgrids,1);
     12md.bed=zeros(md.numberofgrids,1);
     13md.surface=md.thickness;
    2614
    27         disp('      creating temperatures');
    28         md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
     15disp('      creating velocities');
     16md.vx_obs=zeros(md.numberofgrids,1);
     17md.vy_obs=zeros(md.numberofgrids,1);
     18md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    2919
    30         disp('      creating flow law paramter');
    31         md.rheology_B=paterson(md.observed_temperature);
    32         md.rheology_n=3*ones(md.numberofelements,1);
    33        
    34         disp('      creating accumulation rates');
    35         md.accumulation_rate=ones(md.numberofgrids,1)/md.yts; %1m/a
    36         md.melting_rate=0*ones(md.numberofgrids,1)/md.yts; %1m/a
     20disp('      creating drag');
     21md.drag_type=2; %0 none 1 plastic 2 viscous
     22md.drag_coefficient=200*ones(md.numberofgrids,1); %q=1.
     23%Take care of iceshelves: no basal drag
     24pos=find(md.elementoniceshelf);
     25md.drag_coefficient(md.elements(pos,:))=0;
     26md.drag_p=ones(md.numberofelements,1);
     27md.drag_q=ones(md.numberofelements,1);
    3728
    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
     29disp('      creating temperatures');
     30md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
     31
     32disp('      creating flow law paramter');
     33md.rheology_B=paterson(md.observed_temperature);
     34md.rheology_n=3*ones(md.numberofelements,1);
     35
     36disp('      creating accumulation rates');
     37md.accumulation_rate=ones(md.numberofgrids,1)/md.yts; %1m/a
     38md.melting_rate=0*ones(md.numberofgrids,1)/md.yts; %1m/a
     39
     40%Deal with boundary conditions:
     41
     42disp('      boundary conditions for diagnostic model');
     43md=SetMarineIceSheetBC(md,'Front.exp');
     44
     45disp('      boundary conditions for thermal model');
     46md.spctemperature(:,1)=ones(md.numberofgrids,1); %surface temperature
     47md.spctemperature(:,2)=md.observed_temperature;
     48md.geothermalflux=zeros(md.numberofgrids,1);
     49pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2

  • issm/trunk/test/Validation/ThermalTests/Simpleconduction/runme.m

    r2553 r4876  
    2222
    2323%modeled  results
    24 md=solve(md,'analysis_type','thermal','sub_analysis_type','steady');
     24md=solve(md,'analysis_type',ThermalSolutionEnum);
    2525
    2626%plot results
    27 relative=abs((md.results.thermal.temperature-md.observed_temperature)./md.observed_temperature)*100;
    28 relative(find(md.results.thermal.temperature==md.observed_temperature))=0;
    29 plotmodel(md,'data','temperature','title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',md.results.thermal.temperature-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)
     27comp_temp=zeros(md.numberofgrids,1);
     28comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value;
     29relative=abs((comp_temp-md.observed_temperature)./md.observed_temperature)*100;
     30relative(find(comp_temp==md.observed_temperature))=0;
     31plotmodel(md,'data',comp_temp,'title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',comp_temp-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)

  • issm/trunk/test/Validation/ThermalTests/Simpleconduction/thermalboundarycondition.m

    r16 r4876  
    22%Ok, start defining model BC
    33disp('      boundary conditions for thermal model ');
    4 if strcmp(md.type,'3d')
    5         pos1=find(md.elementonbed);     md.gridondirichlet_thermal(md.elements(pos1,1:3))=1; md.dirichletvalues_thermal(md.elements(pos1,1:3))=10;
    6         pos2=find(md.elementonsurface); md.gridondirichlet_thermal(md.elements(pos2,4:6))=1; md.dirichletvalues_thermal(md.elements(pos2,4:6))=0;
     4if md.dim==3,
     5        pos1=find(md.elementonbed);     md.spctemperature(md.elements(pos1,1:3),1)=1; md.spctemperature(md.elements(pos1,1:3),2)=10;
     6        pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6),1)=1; md.spctemperature(md.elements(pos2,4:6),2)=0;
    77        md.vx=zeros(md.numberofgrids,1);
    88        md.vy=zeros(md.numberofgrids,1);

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

    r3768 r4876  
    11%Ok, start defining model parameters here
    22
     3md.dt=0;
    34%parallelization
    4 md.cluster='none';
     5md.np=3;
    56
    67disp('      creating thickness');
    … …  
    3940
    4041disp('      boundary conditions for diagnostic model');
    41 md=SetIceShelfBC(md,'Front.exp');
     42md=SetMarineIceSheetBC(md,'Front.exp');
    4243
    4344disp('      boundary conditions for thermal model');
    44 md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
    45 md.dirichletvalues_thermal=md.observed_temperature;
     45md.spctemperature(:,1)=ones(md.numberofgrids,1); %surface temperature
     46md.spctemperature(:,2)=md.observed_temperature;
    4647md.geothermalflux=zeros(md.numberofgrids,1);
    4748pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=0.1; %100 mW/m^2

  • issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux/runme.m

    r2553 r4876  
    2222
    2323%modeled  results
    24 md=solve(md,'analysis_type','thermal','sub_analysis_type','steady');
     24md=solve(md,'analysis_type',ThermalSolutionEnum);
    2525
    2626%plot results
    27 relative=abs((md.results.thermal.temperature-md.observed_temperature)./md.observed_temperature)*100;
    28 relative(find(md.results.thermal.temperature==md.observed_temperature))=0;
    29 plotmodel(md,'data','temperature','title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',md.results.thermal.temperature-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)
     27comp_temp=zeros(md.numberofgrids,1);
     28comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value;
     29relative=abs((comp_temp-md.observed_temperature)./md.observed_temperature)*100;
     30relative(find(comp_temp==md.observed_temperature))=0;
     31plotmodel(md,'data',comp_temp,'title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',comp_temp-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)

  • issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux/thermalboundarycondition.m

    r16 r4876  
    22%Ok, start defining model BC
    33disp('      boundary conditions for thermal model ');
    4 if strcmp(md.type,'3d')
    5         pos2=find(md.elementonsurface); md.gridondirichlet_thermal(md.elements(pos2,4:6))=1; md.dirichletvalues_thermal(md.elements(pos2,4:6))=0;
     4if md.dim==3,
     5        pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6),1)=1; md.spctemperature(md.elements(pos2,4:6),2)=0;
    66        md.vx=zeros(md.numberofgrids,1);
    77        md.vy=zeros(md.numberofgrids,1);
Note: See TracChangeset for help on using the changeset viewer.