CB_014_xfit_premelt_extended.m

output figures

contents

function VBR = CB_014_xfit_premelt_extended()
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % CB_014_xfit_premelt_extended.m
  %
  %  Calls VBR using xfit_premelt method including the direct melt effect,
  %  from:
  %    Hatsuki Yamauchi and Yasuko Takei, JGR 2024, "Effect of Melt on
  %    Polycrystal Anelasticity", https://doi.org/10.1029/2023JB027738
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  %% write method list %%
  VBR.in.elastic.methods_list={'anharmonic'};
  VBR.in.anelastic.methods_list={'xfit_premelt'};
  VBR.in.anelastic.xfit_premelt.include_direct_melt_effect = 1;

  % load anharmonic parameters, adjust Gu_0_ol and derivatives to match YT2016
  VBR.in.elastic.anharmonic.Gu_0_ol=72.45; %[GPa]
  VBR.in.elastic.anharmonic.dG_dT = -10.94*1e6; % Pa/C    (equivalent ot Pa/K)
  VBR.in.elastic.anharmonic.dG_dP = 1.987; % GPa / GPa

  %% Define the Thermodynamic State %%
  VBR.in.SV.T_K=1200:5:1500;
  VBR.in.SV.T_K=VBR.in.SV.T_K+273;
  sz=size(VBR.in.SV.T_K); % temperature [K]
  VBR.in.SV.P_GPa = full_nd(2.5, sz); % pressure [GPa]

  Tn_cases = [.96, .98, 1.0, 1.1, 1.1, 1.1, 1.1];
  phi_cases = [0., 0., 0., 0.1, 1, 2, 4]*0.01;

  % remaining state variables (ISV)
  VBR.in.SV.dg_um=full_nd(.004*1e6, sz); % grain size [um]
  VBR.in.SV.rho = full_nd(3300, sz); % density [kg m^-3]
  VBR.in.SV.sig_MPa = full_nd(1, sz); % differential stress [MPa]
  VBR.in.SV.f = 1; % 1 Hz

  if ~vbr_tests_are_running()
    f1=figure('PaperPosition',[0,0,6,4],'PaperPositionMode','manual');
    nTn = numel(Tn_cases);    
    for iTn = 1:nTn
          VBRi = VBR;
          VBRi.in.SV.phi = full_nd(phi_cases(iTn), sz);
          VBRi.in.SV.Tsolidus_K = VBR.in.SV.T_K / Tn_cases(iTn);
          [VBRi] = VBR_spine(VBRi) ;

          results.Q = VBRi.out.anelastic.xfit_premelt.Q;          

          dname = [num2str(Tn_cases(iTn)), ', ', num2str(phi_cases(iTn))];
          figure(f1)
          if iTn > 1
              hold all
          end
          plot(VBR.in.SV.T_K - 273, results.Q, 'displayname', dname, 'linewidth', 1.5)
          legend('Location','eastoutside')
          title('Qs vs T for different (Tn, phi) curves')          
          xlabel('Temperature [C]', 'fontsize', 12)
          ylabel('Qs', 'fontsize', 12)
          ylim([0, 200])

          yticks(0:20:200)
    end
    

    saveas(gcf,'./figures/CB_014_xfit_premelt_extended.png')
  end
end