Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Sound velocity anisotropy (Vp, Vs) discrepancy between different versions of MTEX #1830

Closed
ekolesni opened this issue May 17, 2023 · 15 comments
Closed

Sound velocity anisotropy (Vp, Vs) discrepancy between different versions of MTEX #1830

ekolesni opened this issue May 17, 2023 · 15 comments

Comments

@ekolesni
Copy link

What do you want to do?
I would like to calculate sound velocity anisotropy of a powder sample from its IPF
What data do you have?
I use odf obtained from XRD diffraction
What code do you use?
Please provide minimalist code with code in the following form

% load the data
pf = PoleFigure.load(fname,h,CS,SS,'interface','beartex');
odffile_name = calcODF(pf, 'resolution',5*degree)
ipf = Miller(0,0,1,CS)
nextAxis

% extrema Vp
[maxVp_exp, maxVpPos_exp] = max(vp_exp);
[minVp_exp, minVpPos_exp] = min(vp_exp);
% percentage anisotropy Vp
AVp_exp = 200*(maxVp_exp-minVp_exp) / (maxVp_exp+minVp_exp);
nextAxis

% Plot P-wave velocity (km/s)
plotPDF(ebsd,Miller(1,0,0,ebsd.CS)
plot(vp_exp,'contourf','complete','upper')
mtexTitle('Vp (km/s)',titleOpt{:})
xlabel(['Vp Anisotropy = ',num2str(AVp_exp,'%6.1f')],titleOpt{:})

**What result do you get**
I get in 5.9.0 a lower value of Vp anisotropy of powder sample than on earlier version of MTEX (5.3.0)
**What result do you expect**
I am wondering if there is a bug in the last version or there was some problem fixed in previous releases?
**Error Message**
If you get an error message copy and paste it below

put the error message here

**What MTEX version do you use?**
5.9.0
@lfgmorales
Copy link
Contributor

lfgmorales commented May 17, 2023 via email

@ekolesni
Copy link
Author

ekolesni commented May 17, 2023 via email

@Ilya-Ku-I
Copy link

Ilya-Ku-I commented May 30, 2023 •
edited by ralfHielscher

I would like to add several points here. I have tried to trace the origin of the problem and it is not there before version 5.9.0., i.e. before the change of the ODF class to the SO3Fun. Somehow different versions with different classes compute the bulk elastic tensor from ODFs differently. If I use exactly the same input file for ODF and compute the bulk tensor using the same single-crystal Cij tensor I get different results in earlier (before 5.9.0) and later versions.
in earlier versions:

[CVoightFeSi14Textured_exp,CReussFeSi14Textured_exp,CHillFeSi14Textured_exp] = calcTensor(odf141,CFeSi14)
CHillFeSi14Textured_exp = stiffnessTensor (show methods, plot)
  density: 14.597                        
  unit   : GPa                           
  rank   : 4 (3 x 3 x 3 x 3)             
  mineral: Fe (6/mmm, X||a*, Y||b, Z||c*)
 
  tensor in Voigt matrix representation: *10^2
 22.688 12.917  12.13  0.006   0.03  0.058
 12.917 22.562 12.224 -0.006  0.022 -0.003
  12.13 12.224 23.564  0.002 -0.058 -0.044
  0.006 -0.006  0.002  4.832 -0.008  0.003
   0.03  0.022 -0.058 -0.008   4.84  0.014
  0.058 -0.003 -0.044  0.003  0.014  4.983`

in later versions:

[CVoightFeSi14Textured_exp,CReussFeSi14Textured_exp,CHillFeSi14Textured_exp] = calcTensor(odf141,CFeSi14)
CHillFeSi14Textured_exp = stiffnessTensor (xyz)
  density: 14.597           
  unit   : GPa              
  rank   : 4 (3 x 3 x 3 x 3)
 
  tensor in Voigt matrix representation: *10^2
 22.632 12.712 12.419  0.002  0.009  0.021
 12.712 22.582 12.455 -0.002  0.008      0
 12.419 12.455 22.969  0.001  -0.02 -0.016
  0.002 -0.002  0.001  4.968 -0.001      0
  0.009  0.008  -0.02 -0.001  4.976  0.005
  0.021      0 -0.016      0  0.005  4.991

I have tried to play around a bit and here are my observations:
If one saves the odf from earlier versions as a .mat file and then imports it in the later versions (it converts it automatically to SO3Fun class) the computation result in the tensor from the later versions.
If one saves the bulk tensor from earlier versions as a .mat file and then imports it in the later versions the computation result in the tensor from the later versions. Consequently, the velocity calculations are different.

So, once again, for me, it looks like some error from going from the ODF class to the SO3Fun class. The difference in tensors might not look big, but it results in a significant difference in velocity calculations - with a 3-4 times difference in the anisotropy.
I have checked it for several experimental odf-s and the difference is always there. And it is not only a question of which version of mtex is better to use. If the correct calculations are performed only starting from version 5.9.0 it means that all bulk elastic tensors calculated from ODFs using mtex before 2023 are wrong and all publications (including some by our group) using them should be corrected.

@lfgmorales I would greatly appreciate it if you or someone else from the team could have a look at this matter. I am happy to help anyhow if I can.

Best,
Ilya

@lfgmorales
Copy link
Contributor

lfgmorales commented May 30, 2023 via email

@Ilya-Ku-I
Copy link

Hello Luiz,

Thanks for the quick reply! Attached is .zip with odf, single crystal stiffness tensor, bulk stiffness tensor saved in the 5.5.1 version, and the scripts that I use to extract the odf from the experimental data. It is not a ebsd file but a file with Pole Figure Data.

Best,
Ilya
textures-for-tests.zip

@Ilya-Ku-I
Copy link

Ilya-Ku-I commented Jun 29, 2023 •
edited

@lfgmorales
Hello Luiz,

May I ask if you had a chance to have a look?

Best,
Ilya

@lfgmorales
Copy link
Contributor

lfgmorales commented Jun 30, 2023 via email

@kilir
Copy link
Contributor

kilir commented Jul 3, 2023

Hi,
I cannot reproduce your issue. Could you create a reproducible example?
Depending what you do, there have been some bug fixes in the past 3 years, e.g. a4d5c28 related to averaging of densities (of which the issue/discussion I do not appear to be capable to find again here).
Cheers,
Rüdiger

@Ilya-Ku-I
Copy link

@kilir
Hello Rüdiger,

Thanks for the try to solve the issue. Could you please clarify what exactly you cannot reproduce? For me, if I use the same experimental data (or the same EVPSC data - I checked and the mismatch results is the same) to calculate the aggregate tensor I get different results in versions before 5.8.2 (included) or after.

In both versions when I import the data to create the orientation distribution function (odf) and plot the IPF they look identical:
Texture-strong

The odf can be created either from the pole figure data(experimental - via odf = calcODF(pf)) or from Euler angles (EVPSC via odf = calcDensity(data)).
In the earlier versions the resulting file is the 1x1 ODF file while in the later versions the result is the 1x1 SO3FunRBF file.

I use this odf to calculate the bulk tensor. I though that the issue may be related to the high Cij values so I tried with to different single-crystal tensors:
CFeSi14_low = stiffnessTensor (Fe)
density: 9.1
unit : GPa
rank : 4 (3 x 3 x 3 x 3)

tensor in Voigt matrix representation:
438 212 169 0 0 0
212 438 169 0 0 0
169 169 526 0 0 0
0 0 0 100 0 0
0 0 0 0 100 0
0 0 0 0 0 113

CFeSi14 = stiffnessTensor (Fe)
density: 14.597
unit : GPa
rank : 4 (3 x 3 x 3 x 3)

tensor in Voigt matrix representation: *10^2
23.04 13.3 11.29 0 0 0
13.3 23.04 11.29 0 0 0
11.29 11.29 25.52 0 0 0
0 0 0 4.37 0 0
0 0 0 0 4.37 0
0 0 0 0 0 4.87

I do calculations using [CVoightFeSi14Textured,CReussFeSi14Textured,CHillFeSi14Textured] = calcTensor(odf,CFeSi14_low).

The resulting bulk tensors are different in earlier and later versions. In earlier versions with the 1x1 ODF file I get the tensors below and high anisotropy of velocities:

CHillFeSi14Textured = stiffnessTensor (xyz)
density: 9.1
unit : GPa
rank : 4 (3 x 3 x 3 x 3)

tensor in Voigt matrix representation:
433.2 205.7 189.8 0.2 0.6 1.8
205.7 428.5 191.6 -0.1 0.4 0.2
189.8 191.6 464.6 0.1 -1.7 -0.9
0.2 -0.1 0.1 112 -0.2 0.1
0.6 0.4 -1.7 -0.2 112.2 0.4
1.8 0.2 -0.9 0.1 0.4 116
Velocities-Cij-low-5 8 2

CHillFeSi14Textured_exp = stiffnessTensor (xyz)
density: 14.597
unit : GPa
rank : 4 (3 x 3 x 3 x 3)

tensor in Voigt matrix representation: *10^2
22.688 12.917 12.13 0.006 0.03 0.058
12.917 22.562 12.224 -0.006 0.022 -0.003
12.13 12.224 23.564 0.002 -0.058 -0.044
0.006 -0.006 0.002 4.832 -0.008 0.003
0.03 0.022 -0.058 -0.008 4.84 0.014
0.058 -0.003 -0.044 0.003 0.014 4.983
Velocities-Cij-high-5 8 2

In the later versions with the 1x1SO3FunRBF file I get the tensors below and low anisotropy of velocities:
CHillFeSi14Textured = stiffnessTensor (xyz)
density: 9.1
unit : GPa
rank : 4 (3 x 3 x 3 x 3)

tensor in Voigt matrix representation:
433 201.9 196.4 0.1 0.2 0.7
201.9 431 197 0 0.2 0.1
196.4 197 445.4 0.1 -0.6 -0.3
0.1 0 0.1 115.6 0 0
0.2 0.2 -0.6 0 115.8 0.1
0.7 0.1 -0.3 0 0.1 116.2
Velocities-Cij-low-5 10 0

CHillFeSi14Textured_exp = stiffnessTensor (xyz)
density: 14.597
unit : GPa
rank : 4 (3 x 3 x 3 x 3)

tensor in Voigt matrix representation: *10^2
22.633 12.712 12.419 0.002 0.009 0.021
12.712 22.582 12.455 -0.002 0.008 0
12.419 12.455 22.969 0.001 -0.02 -0.016
0.002 -0.002 0.001 4.968 -0.001 0
0.009 0.008 -0.02 -0.001 4.976 0.005
0.021 0 -0.016 0 0.005 4.991
Velocities-Cij-high-5 10 0

These are quite strong textures, so I tried the same with moderate and weak textures - the same problem.
And the topic starter ekolesni has the exact same issue, so the problem is reproducible at different machines.

Moreover, if I save the 1x1 ODF from the earlier versions of mtex as a .mat file and open it in the later version of mtex it automatically converts it to the 1x1SO3FunRBF and the resulting calculations are the same as if I have just created the 1x1SO3FunRBF from the original experimental (or EVPSC) data. I could not check the inverse situation, since earlier versions give an error when trying to open the .mat file with 1x1SO3FunRBF variable.

I am happy to provide files or more examples and would greatly appreciate your help!
Thanks!

Best wishes,
Ilya

@kilir
Copy link
Contributor

kilir commented Jul 13, 2023 •
edited

This one f5bd9f4 introduced the change.

Edit: the odf is fine - I simply forgot to square the norm - is resolved.

vpqtz_odfsynth2_54a5101b1bddbdad35b1a115f0dfad6cdab666e1_version_MTEX 5 8 2_ vpqtz_odfsynth2_f5bd9f470a29471044651b67b83f497ff902f7fe_version_MTEX 5 8 2_ vpqtz_odfsynth2_develop_version_MTEX 5 12 0beta1_ 
cd /path/to/output
commits = { ...
    'cd2db00c57f38c34d7cf2b29643df72838293e93' ...
    '54a5101b1bddbdad35b1a115f0dfad6cdab666e1' ... % large J
    'f5bd9f470a29471044651b67b83f497ff902f7fe' ... % small J <--this one introduced the difference
    '6c86a5fa8be014cce1dd51d2e02374e266ea4b0c' ...
    'develop'
    }
save('commits.mat','commits') % some versions of mtex clear the workspace

for i=1:length(commits)

    save('i.mat','i') % some versions of mtex clear the workspace

    % set desired mtex version
    cd  /path/to/mtexgit
    system(['git checkout ' commits{i}])

    try
        uninstall_mtex
    end

    addpath /path/to/mtexgit
    startup_mtex

    setMTEXpref('figSize','normal');

    fid = fopen('VERSION','r');
    mversion = fgetl(fid);
    fclose(fid);

    % do the tests
    cd /path/to/output
    load('commits.mat')
    load('i.mat')

    mtexdata dubna

    odf = calcODF(pf,'halfwidth',15*degree)
    odf = FourierODF(odf)
    CS = crystalSymmetry('-3m',[  4.9134  4.9134  5.4052],...
        [  90.0000  90.0000 120.0000]*degree,'x||a','z||c');
    % stiffness tensor
    M =....
        [[   86.76    6.868   11.85  -18.02    0.00    0.00];...
        [     6.868   86.76   11.85  +18.02    0.00    0.00];...
        [    11.85    11.85  105.46    0.00    0.00    0.00];...
        [   -18.02   +18.02    0.00   58.14    0.00    0.00];...
        [     0.00     0.00    0.00    0.00   58.14  -18.02];...
        [     0.00     0.00    0.00    0.00  -18.02   39.946]];

    % stiffness tensor
    C = stiffnessTensor(M,CS, 'density' ,2.65);
    C = transformReferenceFrame(C,odf.CS)

    %--------------------------------------------------------------------------

    try
        J = odf.textureindex
    catch
        J = odf.norm ^2
    end

    [CV,CR,CH] = calcTensor(odf,C)
    vp  = CV.velocity;
    anisovp = 200*(max(vp)-min(vp))./(max(vp)+min(vp));

    plot(C,'upper','minmax','complete')
    mtexTitle('single crystal')
    nextAxis
    plotPDF(odf,Miller(0,0,0,1,odf.CS),'minmax','BR',{'J:' num2str(round(J,2))} )
    nextAxis
    plot(vp,'minmax','upper' ,'BR',{'aniso:' num2str(round(anisovp,2))})
    mtexTitle('vp')
    mtexTitle([commits{i} ' ' mversion],'global')
    mtexColorbar

    saveFigure(['vpqtz_odfsynth_' commits{i} '_version_' mversion '_.png'])
end

@ralfHielscher
Copy link
Member

Hi ekolesni,

thank you for reporting this issue. This was indeed a regression against earlier version. Could you please check the fix 6755b39?

All the best,

Ralf.

@lfgmorales
Copy link
Contributor

lfgmorales commented Jul 17, 2023 via email

@ekolesni
Copy link
Author

Hello Ralf,

I checked the updated 6755b39 on 5.10.0 with sqrt, now the calculation result matches with previous versions of MTEX.
Thanks!

Best regards,
Efim

@ralfHielscher
Copy link
Member

By the way, in df09ff3 I was able speed up the computation about a factor of 20. This may not so important for a single average. But if you compute multiple averages it might be helpful.

Ralf.

@ekolesni
Copy link
Author

I tried the df09ff3 on my data, The result is almost the same except a small difference on the order of 10^(-5). Though, for me such level of precision is not important.

Best,
Efim

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
5 participants
@ralfHielscher @kilir @lfgmorales @ekolesni @Ilya-Ku-I