Changeset 25356


Ignore:
Timestamp:
08/07/20 17:56:29 (5 years ago)
Author:
adhikari
Message:

CHG: added a switch that allows us to valideate love_at_depth against analytic solutions for homogeneous case, both for surface loading and tidal forcings. However, got NaN for anlaytic solutions for degrees larger than ~50

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/test/NightlyRun/test2085.m

    r25342 r25356  
    33% Love kernels for degree 2, 20, 200 (tested against analytic solns). 
    44% skip benchmarking for the inner-most interface.
     5
     6% set validation=1 for comparing against the analytic solutions.
     7validation=0;
    58
    69% for volumetric potential
    … …  
    6770        % }}}
    6871
     72% validate tidal potential solutions against the analytic solutions. {{{
     73if validation
     74
     75        param.rho = md.materials.density(1);
     76        param.mu = md.materials.lame_mu(1);
     77        param.G = 6.67e-11;
     78        param.radius = md.materials.radius; 
     79        param.g0 = md.love.g0;
     80        param.source = md.love.forcing_type;
     81
     82        % check against analytic solutions at the following degrees
     83        degrees_analytic = [2 4 8 16 32];
     84        % extract analytic solutions.
     85        for jj=1:length(degrees_analytic)
     86                param.degree = degrees_analytic(jj);
     87                y_temp = yi_analytic_homogenous(param);
     88                if (exist('y_ana','var')~=1)
     89                        y_ana = zeros([length(degrees_analytic) size(y_temp)]);
     90                end
     91                y_ana(jj,:,:)=y_temp;
     92        end
     93        y1_ana=squeeze(y_ana(:,:,1));
     94        y2_ana=squeeze(y_ana(:,:,2));
     95        y3_ana=squeeze(y_ana(:,:,3));
     96        y4_ana=squeeze(y_ana(:,:,4));
     97        y5_ana=squeeze(y_ana(:,:,5));
     98        y6_ana=squeeze(y_ana(:,:,6));
     99
     100        depth = (max(param.radius)-param.radius)/1000; % km.
     101
     102        y1 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1)));
     103        y2 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,2)));
     104        y3 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,3)));
     105        y4 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,4)));
     106        y5 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,5)));
     107        y6 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,6)));
     108       
     109        set(0,'DefaultAxesFontSize',16,'DefaultTextFontSize',15,'DefaultAxesLineWidth',1,...
     110                'DefaultTextInterpreter','Latex','DefaultAxesFontName','Arial','DefaultLineMarkerSize',8)
     111        figure1=figure('Position', [100, 100, 800, 1000]);
     112        %---------------------------------------------------------------------
     113        axes1 = axes('Layer','top','Position',[0.1 0.73 0.4 0.23]);
     114        axis(axes1,[-0.004 0.088 -0.5 6.5]);
     115        ylabel(axes1,'Depth [$\times$1000 km]','fontsize',20);
     116        xlabel(axes1,'$y_1$ [m]','fontsize',20);
     117        box(axes1,'on'); grid(axes1,'on'); hold(axes1,'all');
     118        plot(y1_ana,depth/1000,'linewidth',2,'parent',axes1);
     119        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     120        plot(y1(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes1);
     121        %---------------------------------------------------------------------
     122        axes2 = axes('Layer','top','Position',[0.55 0.73 0.4 0.23]);
     123        axis(axes2,[-200 3200 -0.5 6.5]);
     124        xlabel(axes2,'$y_2$ [Pa]','fontsize',20);
     125        box(axes2,'on'); grid(axes2,'on'); hold(axes2,'all');
     126        plot(y2_ana,depth/1000,'linewidth',2,'parent',axes2);
     127        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     128        plot(y2(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes2);
     129        %---------------------------------------------------------------------
     130        axes3 = axes('Layer','top','Position',[0.1 0.4 0.4 0.23]);
     131        axis(axes3,[-0.002 0.034 -0.5 6.5]);
     132        ylabel(axes3,'Depth [$\times$1000 km]','fontsize',20);
     133        xlabel(axes3,'$y_3$ [m]','fontsize',20);
     134        box(axes3,'on'); grid(axes3,'on'); hold(axes3,'all');
     135        plot(y3_ana,depth/1000,'linewidth',2,'parent',axes3);
     136        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     137        plot(y3(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes3);
     138        %---------------------------------------------------------------------
     139        axes4 = axes('Layer','top','Position',[0.55 0.4 0.4 0.23]);
     140        axis(axes4,[-100 1600 -0.5 6.5]);
     141        xlabel(axes4,'$y_4$ [Pa]','fontsize',20);
     142        box(axes4,'on'); grid(axes4,'on'); hold(axes4,'all');
     143        plot(y4_ana,depth/1000,'linewidth',2,'parent',axes4);
     144        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     145        plot(y4(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes4);
     146        %---------------------------------------------------------------------
     147        axes5 = axes('Layer','top','Position',[0.1 0.07 0.4 0.23]);
     148        axis(axes5,[-0.05 1.5 0 6.5]);
     149        ylabel(axes5,'Depth [$\times$1000 km]','fontsize',20);
     150        xlabel(axes5,'$y_5$ [m$^2$ s$^{-2}$]','fontsize',20);
     151        box(axes5,'on'); grid(axes5,'on'); hold(axes5,'all');
     152        plot(y5_ana,depth/1000,'linewidth',2,'parent',axes5);
     153        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     154        plot(y5(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes5);
     155        %---------------------------------------------------------------------
     156        axes6 = axes('Layer','top','Position',[0.55 0.07 0.4 0.23]);
     157        axis(axes6,[-1e-7 4e-7 0 6.5]);
     158        xlabel(axes6,'$y_6$ [m s$^{-2}$]','fontsize',20);
     159        box(axes6,'on'); grid(axes6,'on'); hold(axes6,'all');
     160        plot(y6_ana,depth/1000,'linewidth',2,'parent',axes6);
     161        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     162        plot(y6(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes6);
     163        %---------------------------------------------------------------------
     164
     165        legend(axes6,'n=2','n=4','n=8','n=16','n=32');
     166        %export_fig('/Users/adhikari/issm/trunk-jpl/src/m/contrib/adhikari/issm_vs_analytic_loading_homogeneous.pdf');
     167       
     168else
     169        %
     170end
     171% }}}
     172
    69173% for surface load.
    70174
    … …  
    97201        y6_loading_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,6)));
    98202        % }}}
     203
     204% validate loading solutions against the analytic solutions. {{{
     205if validation
     206
     207        param.source = md.love.forcing_type;
     208
     209        % extract analytic solutions.
     210        for jj=1:length(degrees_analytic)
     211                param.degree = degrees_analytic(jj);
     212                y_temp = yi_analytic_homogenous(param);
     213                if (exist('y_ana','var')~=1)
     214                        y_ana = zeros([length(degrees_analytic) size(y_temp)]);
     215                end
     216                y_ana(jj,:,:)=y_temp;
     217        end
     218        y1_ana=squeeze(y_ana(:,:,1));
     219        y2_ana=squeeze(y_ana(:,:,2));
     220        y3_ana=squeeze(y_ana(:,:,3));
     221        y4_ana=squeeze(y_ana(:,:,4));
     222        y5_ana=squeeze(y_ana(:,:,5));
     223        y6_ana=squeeze(y_ana(:,:,6));
     224
     225        depth = (max(param.radius)-param.radius)/1000; % km.
     226
     227        y1 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1)));
     228        y2 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,2)));
     229        y3 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,3)));
     230        y4 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,4)));
     231        y5 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,5)));
     232        y6 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,6)));
     233       
     234        set(0,'DefaultAxesFontSize',16,'DefaultTextFontSize',15,'DefaultAxesLineWidth',1,...
     235                'DefaultTextInterpreter','Latex','DefaultAxesFontName','Arial','DefaultLineMarkerSize',8)
     236        figure1=figure('Position', [100, 100, 800, 1000]);
     237        %---------------------------------------------------------------------
     238        axes1 = axes('Layer','top','Position',[0.1 0.73 0.4 0.23]);
     239        axis(axes1,[-0.11 0.01 -0.5 6.5]);
     240        ylabel(axes1,'Depth [$\times$1000 km]','fontsize',20);
     241        xlabel(axes1,'$y_1$ [m]','fontsize',20);
     242        box(axes1,'on'); grid(axes1,'on'); hold(axes1,'all');
     243        plot(y1_ana,depth/1000,'linewidth',2,'parent',axes1);
     244        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     245        plot(y1(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes1);
     246        %---------------------------------------------------------------------
     247        axes2 = axes('Layer','top','Position',[0.55 0.73 0.4 0.23]);
     248        axis(axes2,[-1.8e4 0.2e4 -0.5 6.5]);
     249        xlabel(axes2,'$y_2$ [Pa]','fontsize',20);
     250        box(axes2,'on'); grid(axes2,'on'); hold(axes2,'all');
     251        plot(y2_ana,depth/1000,'linewidth',2,'parent',axes2);
     252        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     253        plot(y2(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes2);
     254        %---------------------------------------------------------------------
     255        axes3 = axes('Layer','top','Position',[0.1 0.4 0.4 0.23]);
     256        axis(axes3,[-0.023 0.002 -0.5 6.5]);
     257        ylabel(axes3,'Depth [$\times$1000 km]','fontsize',20);
     258        xlabel(axes3,'$y_3$ [m]','fontsize',20);
     259        box(axes3,'on'); grid(axes3,'on'); hold(axes3,'all');
     260        plot(y3_ana,depth/1000,'linewidth',2,'parent',axes3);
     261        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     262        plot(y3(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes3);
     263        %---------------------------------------------------------------------
     264        axes4 = axes('Layer','top','Position',[0.55 0.4 0.4 0.23]);
     265        axis(axes4,[-1300 100 -0.5 6.5]);
     266        xlabel(axes4,'$y_4$ [Pa]','fontsize',20);
     267        box(axes4,'on'); grid(axes4,'on'); hold(axes4,'all');
     268        plot(y4_ana,depth/1000,'linewidth',2,'parent',axes4);
     269        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     270        plot(y4(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes4);
     271        %---------------------------------------------------------------------
     272        axes5 = axes('Layer','top','Position',[0.1 0.07 0.4 0.23]);
     273        axis(axes5,[-0.05 0.9 -0.5 6.5]);
     274        ylabel(axes5,'Depth [$\times$1000 km]','fontsize',20);
     275        xlabel(axes5,'$y_5$ [m$^2$ s$^{-2}$]','fontsize',20);
     276        box(axes5,'on'); grid(axes5,'on'); hold(axes5,'all');
     277        plot(y5_ana,depth/1000,'linewidth',2,'parent',axes5);
     278        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     279        plot(y5(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes5);
     280        %---------------------------------------------------------------------
     281        axes6 = axes('Layer','top','Position',[0.55 0.07 0.4 0.23]);
     282        axis(axes6,[-1e-7 1e-6 -0.5 6.5]);
     283        xlabel(axes6,'$y_6$ [m s$^{-2}$]','fontsize',20);
     284        box(axes6,'on'); grid(axes6,'on'); hold(axes6,'all');
     285        plot(y6_ana,depth/1000,'linewidth',2,'parent',axes6);
     286        set(gca,'ColorOrderIndex',1,'Ydir','reverse');
     287        plot(y6(degrees_analytic+1,:),depth/1000,'o','linewidth',2,'parent',axes6);
     288        %---------------------------------------------------------------------
     289
     290        legend(axes6,'n=2','n=4','n=8','n=16','n=32');
     291        %export_fig('/Users/adhikari/issm/trunk-jpl/src/m/contrib/adhikari/issm_vs_analytic_tidal_homogeneous.pdf');
     292       
     293else
     294        %
     295end
     296% }}}
    99297
    100298field_names = {...
Note: See TracChangeset for help on using the changeset viewer.