Skip to content

Commit

Permalink
tests for anisotropy added
Browse files Browse the repository at this point in the history
  • Loading branch information
Bernard Giroux committed Nov 10, 2023
1 parent fd803cd commit 1baddef
Show file tree
Hide file tree
Showing 10 changed files with 445 additions and 61 deletions.
25 changes: 0 additions & 25 deletions tests/files/Grid2Drcdsp_tt_grid_aniso.vtr
View file

This file was deleted.

6 changes: 3 additions & 3 deletions tests/files/Grid2Drcsp_tt_grid_aniso.vtr
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,13 +10,13 @@
<CellData>
</CellData>
<Coordinates>
<DataArray type="Float64" Name="Array 0x11090dfc0" format="binary" RangeMin="0" RangeMax="6000">
<DataArray type="Float64" Name="Array 0x141e04080" format="binary" RangeMin="0" RangeMax="6000">
AQAAAACAAAAoAwAAFwEAAA==eJxNzi1MQmEUxvEbCQYiwUAkGAgG4g0GEiMYSI7AHIndOeacwd0x5phzehVUBJUroCL4cf2EQCAaDEQD4UaDgUgwGJ5/8C2/ved9zzmPYfw/CVOuSHuVe06O1qivS3OT9y3p2/wrSLfI/22ZLtG3I8O79O/Rv4+OHB8w71B6ZeZWpHPE/GNpnbCnKpOn7KvJaJ29ZzJ4zv4LOUW/IbMud0xdkgvjTfJhrEVOjLTJi6ErcmPgmvw4Q+tGfmO6I78weSs/0OzKAUZ7soPhO1nF4L0sofEgN3CK2Ufp44JHHdvo4/yTTGEFxzj3LONYxBH+YuxF5tHDH4y8ygy6OMHQm1xGBz8x8C6X0MYhznCxL62++QcopYc0
</DataArray>
<DataArray type="Float64" Name="Array 0x110904260" format="binary" RangeMin="0" RangeMax="0">
<DataArray type="Float64" Name="Array 0x141e16080" format="binary" RangeMin="0" RangeMax="0">
AQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=
</DataArray>
<DataArray type="Float64" Name="Array 0x1109043a0" format="binary" RangeMin="0" RangeMax="3000">
<DataArray type="Float64" Name="Array 0x141e04300" format="binary" RangeMin="0" RangeMax="3000">
AQAAAACAAACYAQAAqQAAAA==eJxNzikOwkAYhuFKJLISiUQiR6JIJYqgCIo0CASCNIQQQggtW9lpS1ksR+AIHKFHQCIR3ysY82Rm/uWzrP9TNbIuvSb3tnx1eO9K0+O/LzOPuoGMh9SPZGNM30QWpvTP6PcxkO858xbyuWTuSgZr5ofS3bBnK50d+/aydGDvUeZP7D/LD2aRbMXcsZaQCysX8mE5JScWr+RF+0ZuzN3Jj190H+YHwvo/dQ==
</DataArray>
</Coordinates>
Expand Down
38 changes: 38 additions & 0 deletions tests/files/Grid2Ducsp_tt_grid_aniso.vtu
View file

Large diffs are not rendered by default.

39 changes: 23 additions & 16 deletions tests/files/elliptical_fine2d.vtr
View file
Original file line number Diff line number Diff line change
@@ -1,21 +1,28 @@
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
<RectilinearGrid WholeExtent="0 100 0 0 0 50">
<Piece Extent="0 100 0 0 0 50" >
<PointData>
</PointData>
<CellData>
<DataArray type="Float64" Name="Slowness" format="appended" RangeMin="0.00016184936638" RangeMax="0.000344114832" offset="0" />
<DataArray type="Float64" Name="xi" format="appended" RangeMin="1.2649110641" RangeMax="1.2649110641" offset="820" />
</CellData>
<Coordinates>
<DataArray type="Float32" Name="Array 0x14422f430" format="appended" RangeMin="" RangeMax="" offset="1132" />
<DataArray type="Float32" Name="Array 0x144229fc0" format="appended" RangeMin="" RangeMax="" offset="1500" />
<DataArray type="Float32" Name="Array 0x144212ba0" format="appended" RangeMin="" RangeMax="" offset="1540" />
</Coordinates>
</Piece>
<Piece Extent="0 100 0 0 0 50">
<PointData>
</PointData>
<CellData>
<DataArray type="Float64" Name="Slowness" format="binary" RangeMin="0.00016184936638457496" RangeMax="0.0003441148320023047">
AgAAAACAAABAHAAAwgEAAI4AAAA=eJztxU1L0wEAx/ERNLCJFV5knUIsXTb17zZhMg0SQ71YBwtRzMQRjdT5dBBKDQyMwoiCSYdA8aHyofmY5xAWJbLUGMOyEs1LGE5PIvoC+r2Dvp/LJ9KfGCt64sqPMPM/d8X9/fzCwszqaxvexec9TmYW749OrW4lMbP6wP/4Q+6Qg5nFr8eO1+W5mVn9zfiRNhzOZmbxV+Ol2XWXmdWPSm4XvolnZvWztda1T+MGM4ud9c0Wbxkzqy2/6tavmJhZHTN3mwtGs5hZ/D3FVpxZycxqq7uzceMkM6u3Fit28uYzmVnsCrUFrfeYWe1ufhq5kcPM6oSPoYLZ3xnMLE51hm8l1DKzui975MuDTTszi0ds04F1HzOrrWf+zBzbu8jM4oc/fbGpTmZWl5u2+8OnmVltP2e6bwykM7M4+YTHN5nLzOp3g9erjZULzCzu7vV7evzMrPa8b+97e4qZ1UvbvuidCRszi0taJueCZcysTi8u7PIepDGzuOZyU7RhkJnVdkcwMHeVmdXtE9YF22EqM4t3s1pfjY8xs3r5cOiS4yYzq4vi66MDicysrmoo7dgJnWdm8XBp4GxcB/P/+RGVow7DeJxbETBDkbNB3X4FDtpAObLmpsUoPUqP0tjoh60hgiGf1UbpUXqUxkJz798lkbV+lB6lR2lsdP+e9eW8uaP0KD1KY6PZP7duM9AZpUfpURobfe/xgu/H3qiO0qP0KI2FFg/+9O7wulF6lB6lsdHLw9luSReN0qP0KI2NViifXHHYbJQepUdpbDQApqAU2Q==
</DataArray>
<DataArray type="Float64" Name="xi" format="binary" RangeMin="1.2649110640673515" RangeMax="1.264911064067352">
AgAAAACAAABAHAAAlAAAAD8AAAA=eJzt2rENgEAMA8ANeCTmQbSwBvuXDMAEIEHjD1x1E1iOkrRt3af5WBpJkg8cSV46kCT5wvR8R/ZsOp8kyZqm+4skSZLkf0zfJ8meTeeTJFnT9H6ZJFnTdH+RJEmSJEmSJN+Z/m8hSdY03V8kSZIkSZLknen/epJkTdP9RZKsabq/yJ5N55MkSfJrpuc7MukJNGFQBXic7dcxDQAgDABBB0CCHsIKNvA/IgAFnRuSm07BD1/2On3cWQIrSZIkSZJMN/p2knc2kqHZfZIk//QBQpzlVQ==
</DataArray>
</CellData>
<Coordinates>
<DataArray type="Float64" Name="Array 0x1306073f0" format="binary" RangeMin="0" RangeMax="6000">
AQAAAACAAAAoAwAAFwEAAA==eJxNzi1MQmEUxvEbCQYiwUAkGAgG4g0GEiMYSI7AHIndOeacwd0x5phzehVUBJUroCL4cf2EQCAaDEQD4UaDgUgwGJ5/8C2/ved9zzmPYfw/CVOuSHuVe06O1qivS3OT9y3p2/wrSLfI/22ZLtG3I8O79O/Rv4+OHB8w71B6ZeZWpHPE/GNpnbCnKpOn7KvJaJ29ZzJ4zv4LOUW/IbMud0xdkgvjTfJhrEVOjLTJi6ErcmPgmvw4Q+tGfmO6I78weSs/0OzKAUZ7soPhO1nF4L0sofEgN3CK2Ufp44JHHdvo4/yTTGEFxzj3LONYxBH+YuxF5tHDH4y8ygy6OMHQm1xGBz8x8C6X0MYhznCxL62++QcopYc0
</DataArray>
<DataArray type="Float64" Name="Array 0x130607530" format="binary" RangeMin="0" RangeMax="0">
AQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=
</DataArray>
<DataArray type="Float64" Name="Array 0x130608990" format="binary" RangeMin="0" RangeMax="3000">
AQAAAACAAACYAQAAqQAAAA==eJxNzikOwkAYhuFKJLISiUQiR6JIJYqgCIo0CASCNIQQQggtW9lpS1ksR+AIHKFHQCIR3ysY82Rm/uWzrP9TNbIuvSb3tnx1eO9K0+O/LzOPuoGMh9SPZGNM30QWpvTP6PcxkO858xbyuWTuSgZr5ofS3bBnK50d+/aydGDvUeZP7D/LD2aRbMXcsZaQCysX8mE5JScWr+RF+0ZuzN3Jj190H+YHwvo/dQ==
</DataArray>
</Coordinates>
</Piece>
</RectilinearGrid>
<AppendedData encoding="base64">
_AgAAAACAAABAHAAAwgEAAI4AAAA=eJztxU1L0wEAx/ERNLCJFV5knUIsXTb17zZhMg0SQ71YBwtRzMQRjdT5dBBKDQyMwoiCSYdA8aHyofmY5xAWJbLUGMOyEs1LGE5PIvoC+r2Dvp/LJ9KfGCt64sqPMPM/d8X9/fzCwszqaxvexec9TmYW749OrW4lMbP6wP/4Q+6Qg5nFr8eO1+W5mVn9zfiRNhzOZmbxV+Ol2XWXmdWPSm4XvolnZvWztda1T+MGM4ud9c0Wbxkzqy2/6tavmJhZHTN3mwtGs5hZ/D3FVpxZycxqq7uzceMkM6u3Fit28uYzmVnsCrUFrfeYWe1ufhq5kcPM6oSPoYLZ3xnMLE51hm8l1DKzui975MuDTTszi0ds04F1HzOrrWf+zBzbu8jM4oc/fbGpTmZWl5u2+8OnmVltP2e6bwykM7M4+YTHN5nLzOp3g9erjZULzCzu7vV7evzMrPa8b+97e4qZ1UvbvuidCRszi0taJueCZcysTi8u7PIepDGzuOZyU7RhkJnVdkcwMHeVmdXtE9YF22EqM4t3s1pfjY8xs3r5cOiS4yYzq4vi66MDicysrmoo7dgJnWdm8XBp4GxcB/P/+RGVow7DeJxbETBDkbNB3X4FDtpAObLmpsUoPUqP0tjoh60hgiGf1UbpUXqUxkJz798lkbV+lB6lR2lsdP+e9eW8uaP0KD1KY6PZP7duM9AZpUfpURobfe/xgu/H3qiO0qP0KI2FFg/+9O7wulF6lB6lsdHLw9luSReN0qP0KI2NViifXHHYbJQepUdpbDQApqAU2Q==AgAAAACAAABAHAAAlAAAAD8AAAA=eJzt2rENgEAMA8ANeCTmQbSwBvuXDMAEIEHjD1x1E1iOkrRt3af5WBpJkg8cSV46kCT5wvR8R/ZsOp8kyZqm+4skSZLkf0zfJ8meTeeTJFnT9H6ZJFnTdH+RJEmSJEmSJN+Z/m8hSdY03V8kSZIkSZLknen/epJkTdP9RZKsabq/yJ5N55MkSfJrpuc7MukJNGFQBXic7dcxDQAgDABBB0CCHsIKNvA/IgAFnRuSm07BD1/2On3cWQIrSZIkSZJMN/p2knc2kqHZfZIk//QBQpzlVQ==AQAAAACAAACUAQAAAgEAAA==eJwVzTEoBFAcx/F/GZTBYLgyKIYrg8FwZTC43D2DwXBlMChXBoPhymAwIOkkXZJOkpBO0iXpJN2ibMbbbGzG24x8vPr1fd++r17E/6lMR3QtV3C3urWsY13rLUZkLGs5K1jJylaxTasWY6uGdTzHBjaxhW18w3fs4Ad+4Td28Qd/i5HvSf5K8dqXYqs/RX6AZ/ggH+IjPMtH+Rgf5zk+wSf5FC/wGT7L53iJz/MFvsjLfJmv8FVe4Wt8nW/wzRRL2+47KYZ3eTXF5543+ykuD7yr6Yf6kX6s1/UT/VQ/08/1C/1Kv9Yb+o1+q9/pTf1ef9Af9Zb+pD/rL3o7/QFQUmm9AQAAAACAAAAEAAAADAAAAA==eJxjYGBgAAAABAABAQAAAACAAADMAAAAkgAAAA==eJwVzaEOglAYhuFvM7AZDAY2g8HAZjAYzmagMDn/DXgH3gF3II1GpBlpRhrFjUakGW3G04z6crZ3z/+lI82vOEuBXM5NDXU0UaDISzEl5CinC12poBtVXmWNDd6xxQd22OOAI074wjd+MOAXf17ZwvjL9FyaypUpW7Nj9oa9Ze/YCXvPPrCPbMc+sVP7A0DFLZk=
</AppendedData>
</VTKFile>
41 changes: 41 additions & 0 deletions tests/files/elliptical_fine2d.vtu
View file

Large diffs are not rendered by default.

241 changes: 241 additions & 0 deletions tests/files/prepare_tests_aniso.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,241 @@
#! /usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import vtk


# %%

x = np.arange(0.0, 6000.01, 60.0)
z = np.arange(0.0, 3000.01, 60.0)


a = 2271.0
b = 0.88
chi = 0.3

with open('rcv2daniso.dat', 'w') as f:
print('{0:d}'.format(x.size*z.size), file=f)
for ix in x:
for iz in z:
print('{0:g} {1:g}'.format(ix, iz), file=f)


# %%

x[0] = 0.001

xx, zz = np.meshgrid(x, z)

p = 2*xx / np.sqrt((xx*xx + (1+2*chi) * zz*zz) *
((2*a + b*zz)**2 * (1+2*chi) + b*b * xx*xx))

# p[0, 0] = 0.0

t = 1/b * (np.arctanh(p * b * xx - np.sqrt(1 - (1+2*chi) * p*p * a*a)) +
np.arctanh(np.sqrt(1 - (1+2*chi) * p*p * a*a)))

# t[0, 0] = 0
# t[:, 0] = z / vz


# %%

h = plt.pcolor(xx, zz, t)
plt.contour(xx, zz, t, colors='k', linewidths=0.5)
plt.gca().invert_yaxis()
plt.gca().axis('equal')
plt.colorbar(h)
plt.tight_layout()
plt.show()

# %%

x = np.arange(30.0, 6000.0, 60.0)
z = np.arange(30.0, 3000.0, 60.0)

xx, zz = np.meshgrid(x, z)

vz = a + b*zz
vx = np.sqrt( chi * 2 * vz**2 + vz**2)


sx = 1/vx
sz = 1/vz

xi = sz / sx

plt.imshow(xi)
plt.colorbar()
plt.show()

plt.imshow(sx)
plt.colorbar()
plt.show()



xCoords = vtk.vtkDoubleArray()
yCoords = vtk.vtkDoubleArray()
zCoords = vtk.vtkDoubleArray()

for i in np.arange(0.0, 6000.1, 60.0):
xCoords.InsertNextValue(i)
for i in np.arange(0.0, 3000.1, 60.0):
zCoords.InsertNextValue(i)
yCoords.InsertNextValue(0.0)

rgrid = vtk.vtkRectilinearGrid()
rgrid.SetDimensions(xCoords.GetNumberOfValues(),
yCoords.GetNumberOfValues(),
zCoords.GetNumberOfValues())
rgrid.SetXCoordinates(xCoords)
rgrid.SetYCoordinates(yCoords)
rgrid.SetZCoordinates(zCoords)

slown = vtk.vtkDoubleArray()
for i in sx.flatten():
slown.InsertNextValue(i)
slown.SetName("Slowness")
rgrid.GetCellData().AddArray(slown)

xi_vtk = vtk.vtkDoubleArray()
for i in xi.flatten():
xi_vtk.InsertNextValue(i)
xi_vtk.SetName("xi")
rgrid.GetCellData().AddArray(xi_vtk)

writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetInputData(rgrid)
writer.SetFileName('elliptical_fine2d.vtr')
writer.SetDataModeToBinary()
writer.Write()

# %%
rgrid = vtk.vtkRectilinearGrid()
rgrid.SetDimensions(xCoords.GetNumberOfValues(),
yCoords.GetNumberOfValues(),
zCoords.GetNumberOfValues())
rgrid.SetXCoordinates(xCoords)
rgrid.SetYCoordinates(yCoords)
rgrid.SetZCoordinates(zCoords)

tt = vtk.vtkDoubleArray()
for i in t.flatten():
tt.InsertNextValue(i)
tt.SetName("Travel Time")
rgrid.GetPointData().AddArray(tt)

writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetInputData(rgrid)
writer.SetFileName('sol_analytique_elliptical_2d_tt.vtr')
writer.SetDataModeToBinary()
writer.Write()

# %%

x = np.arange(0.0, 6000.01, 60.0)
z = np.arange(0.0, 3000.01, 60.0)

xx, zz = np.meshgrid(x, z)

xx = xx.flatten()
zz = zz.flatten()
t = t.flatten()

npx = x.size
ncx = x.size - 1
ncz = z.size - 1

tri = []
slo = []
xi = []
for nz in range(ncz):
for nx in range(ncx):
tri.append([nz*npx + nx, nz*npx + nx+1, (nz+1)*npx + nx])

zc = (zz[tri[-1][0]] + zz[tri[-1][1]] + zz[tri[-1][2]]) / 3.
vz = a + b*zc
vx = np.sqrt( chi * 2 * vz**2 + vz**2)

sx = 1/vx
sz = 1/vz

slo.append(sx)
xi.append( sz / sx )

tri.append([nz*npx + nx+1, (nz+1)*npx + nx, (nz+1)*npx + nx + 1])

zc = (zz[tri[-1][0]] + zz[tri[-1][1]] + zz[tri[-1][2]]) / 3.
vz = a + b*zc
vx = np.sqrt( chi * 2 * vz**2 + vz**2)

sx = 1/vx
sz = 1/vz

slo.append(sx)
xi.append( sz / sx )

# %%

ugrid = vtk.vtkUnstructuredGrid()
pts = vtk.vtkPoints()

for n in range(len(xx)):
pts.InsertNextPoint(xx[n], 0.0, zz[n])

ugrid.SetPoints(pts)

triangle = vtk.vtkTriangle()
d_s = vtk.vtkDoubleArray()
d_x = vtk.vtkDoubleArray()
d_s.SetName("Slowness")
d_x.SetName("xi")
for n in range(len(tri)):
triangle.GetPointIds().SetId(0, tri[n][0])
triangle.GetPointIds().SetId(1, tri[n][1])
triangle.GetPointIds().SetId(2, tri[n][2])
ugrid.InsertNextCell( triangle.GetCellType(), triangle.GetPointIds() )
d_s.InsertNextValue(slo[n])
d_x.InsertNextValue(xi[n])

ugrid.GetCellData().AddArray(d_s)
ugrid.GetCellData().AddArray(d_x)

writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetInputData(ugrid)
writer.SetFileName('elliptical_fine2d.vtu')
writer.SetDataModeToBinary()
writer.Write()

# %%

ugrid = vtk.vtkUnstructuredGrid()
pts = vtk.vtkPoints()

for n in range(len(xx)):
pts.InsertNextPoint(xx[n], 0.0, zz[n])

ugrid.SetPoints(pts)

triangle = vtk.vtkTriangle()
for n in range(len(tri)):
triangle.GetPointIds().SetId(0, tri[n][0])
triangle.GetPointIds().SetId(1, tri[n][1])
triangle.GetPointIds().SetId(2, tri[n][2])
ugrid.InsertNextCell( triangle.GetCellType(), triangle.GetPointIds() )

d_t = vtk.vtkDoubleArray()
d_t.SetName("Travel Time")
for n in range(len(t)):
d_t.InsertNextValue(t[n])

ugrid.GetPointData().AddArray(d_t)

# writer = vtk.vtkXMLUnstructuredGridWriter()
# writer.SetInputData(ugrid)
# writer.SetFileName('tests/files/sol_analytique_elliptical_2d_tt.vtu')
# writer.SetDataModeToBinary()
# writer.Write()
35 changes: 20 additions & 15 deletions tests/files/sol_analytique_elliptical_2d_tt.vtr
View file

Large diffs are not rendered by default.

6 changes: 4 additions & 2 deletions tests/test_grid2d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,7 +151,8 @@ const char* models_L[] = {
"./files/layers_coarse2d.vtu"
};
const char* models_aniso[] = {
"./files/elliptical_fine2d.vtr" //,
"./files/elliptical_fine2d.vtr",
"./files/elliptical_fine2d.vtu" //,
// "./files/weakly_an_fine2d.vtu"
};
const char* references[] = {
Expand All @@ -165,6 +166,7 @@ const char* references_L[] = {
"./files/sol_analytique_couches2d_tt.vtr"
};
const char* references_aniso[] = {
"./files/sol_analytique_elliptical_2d_tt.vtr",
"./files/sol_analytique_elliptical_2d_tt.vtr"//,
// "./files/sol_analytique_weakly_an_2d_tt.vtr"
};
Expand Down Expand Up @@ -380,6 +382,6 @@ BOOST_DATA_TEST_CASE(
double error = get_rel_error(ref, rcv);
BOOST_TEST_MESSAGE( "\t\t" << get_class_name(g) << " - error = " << error );

BOOST_TEST(error < 0.01);
BOOST_TEST(error < 0.005);
}

33 changes: 33 additions & 0 deletions tests/test_rgrid2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,39 @@ def test_Grid2Ddsp(self):
self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01,
'DSPM accuracy failed (slowness at nodes)')

class TestAniso(unittest.TestCase):

def setUp(self):
reader = vtk.vtkXMLRectilinearGridReader()
reader.SetFileName('./files/elliptical_fine2d.vtr')
reader.Update()

data = reader.GetOutput()
self.x = vtk_to_numpy(data.GetXCoordinates())
self.z = vtk_to_numpy(data.GetZCoordinates())

self.slowness = vtk_to_numpy(data.GetCellData().GetArray('Slowness'))
dim = (self.x.size-1, self.z.size-1)
self.slowness = self.slowness.reshape(dim, order='F').flatten()

self.xi = vtk_to_numpy(data.GetCellData().GetArray('xi'))
self.xi = self.xi.reshape(dim, order='F').flatten()

self.src = np.loadtxt('./files/src2d.dat', skiprows=1)
self.src = np.roll(self.src, 1).reshape((1, 3))
self.rcv = np.loadtxt('./files/rcv2daniso.dat', skiprows=1)

def test_Grid2Dsp(self):
g = rg.Grid2d(self.x, self.z, method='SPM', nsnx=10, nsnz=10, aniso='elliptical')
g.set_slowness(self.slowness)
g.set_xi(self.xi)
tt = g.raytrace(self.src, self.rcv)
tt = g.get_grid_traveltimes()
tt = tt.flatten()
tt_ref = get_tt('./files/Grid2Drcsp_tt_grid_aniso.vtr')
self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01,
'SPM accuracy failed (elliptical anisotropy)')


class Data_kernel(unittest.TestCase):

Expand Down
Loading

0 comments on commit 1baddef

Please sign in to comment.