Changeset 24775


Ignore:
Timestamp:
05/04/20 09:56:20 (5 years ago)
Author:
Eric.Larour
Message:

CHG: slr test out of a sealevel class model, concept is working.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/test/SandBox/test2004.m

    r24764 r24775  
    77%}}}
    88
    9 skipmesh=1;
    10 skipparameterization=0;
    11 
    129%create sealevel model to hold our information:
    13 if ~skipmesh,
    1410        sl=sealevelmodel();
    15 %Create basins using boundaries from shapefile: %{{{
     11
     12%Create basins using boundaries from shapefile:
     13%some projections we'll rely on:  %{{{
     14proj4326=epsg2proj(4326);
     15proj3031=epsg2proj(3031);
     16%}}}
    1617%HemisphereWest: {{{
    17         sl.addbasin(basin('continent','hemispherewest','name','hemispherewest','epsg',3587,'boundaries',{... %Peru projection 3587
    18         boundary('shppath',shppath,'shpfilename','HemisphereSplit','epsg',4326,'orientation','reverse'),...
    19         boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031),...
    20         boundary('shppath',shppath,'shpfilename','RonneBrunt','epsg',3031,'orientation','reverse'),...
    21         boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031),...
    22         boundary('shppath',shppath,'shpfilename','RonneFront','epsg',3031,'orientation','reverse'),...
    23         boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031),...
    24         boundary('shppath',shppath,'shpfilename','WestAntarctica2','epsg',3031,'orientation','reverse'),...
    25         boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031)...
     18sl.addbasin(basin('continent','hemispherewest','name','hemispherewest','proj',laea(0,-90),'boundaries',{... %Peru projection 3587
     19        boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326,'orientation','reverse'),...
     20        boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),...
     21        boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031,'orientation','reverse'),...
     22        boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),...
     23        boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse'),...
     24        boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),...
     25        boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031,'orientation','reverse'),...
     26        boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031)...
    2627        }));
    2728        %}}}
    2829        %HemisphereEast: {{{
    29         sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','epsg',4462,'boundaries',{... %Australian projection lat,long
    30         boundary('shppath',shppath,'shpfilename','HemisphereSplit','epsg',4326),...
    31         boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031),...
    32         boundary('shppath',shppath,'shpfilename','EastAntarctica2','epsg',3031,'orientation','reverse'),...
    33         boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031)...
     30        sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long
     31                boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),...
     32                boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
     33                boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),...
     34                boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)...
     35                }));
     36%}}}
     37%Antarctica excluding Ronne: {{{
     38sl.addbasin(basin('continent','antarctica','name','antarctica-grounded','proj',proj3031,'boundaries',{...
     39        boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),...
     40        boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031),...
     41        boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
     42        boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031)...
     43        boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031)...
     44        boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031)...
     45        boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031)...
     46        boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031)...
    3447        }));
    3548        %}}}
    36         %Antarctica excluding Ronne: {{{
    37         sl.addbasin(basin('continent','antarctica','name','antarctica-grounded','epsg',3031,'boundaries',{...
    38         boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031),...
    39         boundary('shppath',shppath,'shpfilename','EastAntarctica2','epsg',3031),...
    40         boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031),...
    41         boundary('shppath',shppath,'shpfilename','WestAntarctica2','epsg',3031)...
    42         boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031)...
    43         boundary('shppath',shppath,'shpfilename','RonneIceShelf','epsg',3031)...
    44         boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031)...
    45         boundary('shppath',shppath,'shpfilename','RonneBrunt','epsg',3031)...
    46         }));
    47         %}}}
    4849        %Ronne: {{{
    49         sl.addbasin(basin('continent','antarctica','name','ronne','epsg',3031,'boundaries',{...
    50         boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031),...
    51         boundary('shppath',shppath,'shpfilename','RonneIceShelf','epsg',3031),...
    52         boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031),...
    53         boundary('shppath',shppath,'shpfilename','RonneFront','epsg',3031,'orientation','reverse')...
    54         }));
    55         %}}}
    56 %}}}
     50        sl.addbasin(basin('continent','antarctica','name','ronne','proj',proj3031,'boundaries',{...
     51                boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),...
     52                boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),...
     53                boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),...
     54                boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse')...
     55                }));
     56%}}}
     57
     58%Meshing
    5759%Go through basins and mesh:  %{{{
    5860%meshing parameters:  {{{
    59 hmin=500; hmax=2000; hmin=hmin*1000; hmax=hmax*1000;
     61hmin=2000; hmax=10000; hmin=hmin*1000; hmax=hmax*1000;
    6062tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes
    61 threshold=1;
     63threshold=5;
    6264defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinment',1};
    6365alreadyloaded=0;
    … …  
    7981
    8082        %miscellaneous:
    81         md.mesh.epsg=bas.epsg; md.miscellaneous.name=bas.name;
     83        md.mesh.proj=bas.proj; md.miscellaneous.name=bas.name;
    8284
    8385        %recover mask where we have land:
    … …  
    9193end
    9294%}}}
    93 end
    94 if ~skipparameterization,
     95
     96%Parameterization:
    9597%Parameterize ice sheets : {{{
    9698
    … …  
    123125        end
    124126        if bas.isnameany('ronne'),
    125                
     127
    126128                md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
    127129
    … …  
    133135        %}}}
    134136        %latlong:  % {{{
    135         [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326');
     137        [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326');
    136138        %}}}
    137139        %geometry: {{{
    … …  
    190192
    191193        %recover lat,long:
    192         [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326');
     194        [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326');
    193195
    194196        %mask:  %{{{
    … …  
    244246        %grounded ice:
    245247        md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
    246        
     248
    247249        md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
    248250
    … …  
    252254
    253255        md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
     256        md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
    254257        md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
    255258        md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
    … …  
    275278end
    276279% }}}
    277 end
    278 %Assemble Earth in 3D {{{
     280
     281%%Assemble Earth in 3D {{{
    279282
    280283%parameters:
    … …  
    287290
    288291%figure out how each icecap's mesh connects to the larger earth mesh:
    289 sl.intersections();
     292sl.intersections('force',1);
    290293
    291294%figure out connectivity:
    … …  
    307310sl.transfer('mesh.long');
    308311sl.transfer('slr.deltathickness');
     312sl.transfer('slr.spcthickness');
    309313sl.transfer('slr.Ngia');
    310314sl.transfer('slr.Ugia');
    311 error;
     315sl.transfer('slr.hydro_rate');
     316sl.transfer('slr.sealevel');
     317sl.transfer('dsl.sea_surface_height_change_above_geoid');
     318sl.transfer('dsl.sea_water_pressure_change_at_sea_floor');
    312319
    313320%radius:
    … …  
    315322
    316323%check on the mesh transitions: {{{
     324plotting=1;
    317325if plotting,
    318         flags=ones(sl.earth.mesh.numberofvertices,1);
    319         for i=1:length(sl.transitions),
    320                 flags(sl.transitions{i})=i;
     326        flags=ones(sl.earth.mesh.numberofelements,1);
     327        for i=1:length(sl.eltransitions),
     328                flags(sl.eltransitions{i})=i;
    321329        end
    322330        plotmodel(sl.earth,'data',flags,'shading','faceted','coastline','on','coast_color','g')
    … …  
    325333
    326334% }}}
    327 error;
     335
    328336%Solve Sea-level equation on Earth only:  {{{
    329337md=sl.earth; %we don't do computations on ice sheets or land.
     338
     339%Materials:
     340md.materials=materials('hydro');
    330341
    331342%elastic loading from love numbers:
    … …  
    340351
    341352%New stuff
    342 md.slr.spcthickness = NaN(md.mesh.numberofvertices,1);
    343 md.slr.Ngia = zeros(md.mesh.numberofvertices,1);
    344 md.slr.Ugia = zeros(md.mesh.numberofvertices,1);
     353md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    345354
    346355%Solution parameters
    … …  
    349358md.slr.geodetic=1;
    350359
     360%eustatic:
     361md.slr.rigid=0; md.slr.elastic=0; md.slr.rotation=0;
     362md.cluster=generic('name',oshostname(),'np',10);
     363md=solve(md,'Sealevelrise');
     364
    351365%eustatic + rigid + elastic run:
    352366md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=0;
    353 md.cluster=generic('name',oshostname(),'np',3);
    354 %md.verbose=verbose('111111111');
     367md.cluster=generic('name',oshostname(),'np',16);
     368md.verbose=verbose('111111111');
    355369md=solve(md,'Sealevelrise');
    356370SnoRotation=md.results.SealevelriseSolution.Sealevel;
Note: See TracChangeset for help on using the changeset viewer.