Skip to content

Commit

Permalink
Merge pull request #62 from numericalEFT/compositemesh
Browse files Browse the repository at this point in the history 
Bug fix.
  • Loading branch information
@iintSjds
iintSjds committed Nov 29, 2022
2 parents f9a04d6 + 2c888f2 commit 27b5227
Show file tree
Hide file tree
Showing 4 changed files with 91 additions and 134 deletions.
173 changes: 42 additions & 131 deletions example/compositemesh.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,146 +7,57 @@ using BrillouinZoneMeshes.AbstractMeshes
using BrillouinZoneMeshes.LinearAlgebra
using BrillouinZoneMeshes.StaticArrays

using Test
lattice = Matrix([1 0; 0 1.0]')
atoms = [1]
positions = [zeros(2)]
br = BZMeshes.Cell(lattice=lattice, atoms=atoms, positions=positions)

_dispersion(k) = -sum(cos.(k)) + 0.1
# _dispersion(k) = dot(k, k) - π^2 / 4
T = 0.001

function dispersion(k)
# constrain to 1st bz
_k = [(abs(p) < π) ? p : π for p in k]
return _dispersion(_k)
end

function AbstractMeshes.interval(mesh::ProdMesh, I::Int)
inds = AbstractMeshes._ind2inds(size(mesh), I)
J = Base._sub2ind(size(mesh)[2:end], inds[2:end]...)
x1, x2 = AbstractMeshes.interval(mesh.grids[J], inds[1])
if length(inds) > 2
return [(x1, x2), AbstractMeshes.interval(mesh.mesh, J)...]
function integrand(k)
#if dot(k, k) <= π^2
if abs(k[1]) <= π && abs(k[2]) <= π
return 1 / (dispersion(k)^2 + T^2)
else
return [(x1, x2), AbstractMeshes.interval(mesh.mesh, J)]
return 0.0
end
end

struct CompositeMesh{T,PM,SM} <: AbstractMesh{T,1}
panelmesh::PM
submeshes::Vector{SM}
end

function CompositeMesh(panelmesh::PM, N) where {PM}
T = eltype(panelmesh)
submeshes = []
for (i, p) in enumerate(panelmesh)
intervals = AbstractMeshes.interval(panelmesh, i)
DIM = length(intervals)
origin = [intervals[j][1] for j in 1:DIM]
lattice = diagm(DIM, DIM, [intervals[j][2] - intervals[j][1] for j in 1:DIM])
cm = ChebMesh(origin, lattice, DIM, N)
push!(submeshes, cm)
function int_mesh(mesh)
data = zeros(Float64, size(mesh))
for i in 1:length(mesh)
data[i] = integrand(mesh[i])
end

SM = typeof(submeshes[1])
return CompositeMesh{T,PM,SM}(panelmesh, submeshes)
return integrate(data, mesh)
end

N = 8
bound = [-π, π]
theta = SimpleGrid.Uniform(bound, N; isperiodic=true)
# bzmesh = BZMeshes.PolarMesh(dispersion=dispersion, anglemesh=theta, cell=br,
# kmax=π, Nloggrid=5, Nbasegrid=4, minterval=0.001)
pm = BZMeshes.CompositePolarMesh(dispersion=dispersion,
anglemesh=theta, cell=br, basegridtype=:uniform,
kmax=π * sqrt(2.4), Nloggrid=8, Nbasegrid=3, minterval=0.1T, N=4)

@testset "CompositeMesh" begin
@testset "Constructor" begin
a, b = 0.8, 1.2

N, M = 3, 2
# theta grid dense around 0 and π
theta = CompositeGrid.LogDensedGrid(
:cheb,
[-π, π],
[-π, 0, π],
N,
0.1,
M
)
println(theta)
grids = [CompositeGrid.LogDensedGrid(:cheb, [0.0, 2.0], [sqrt(a * cos(θ)^2 + b * sin(θ)^2),], N, 0.1, M) for θ in theta]

pm = ProdMesh(theta, grids)
cm = CompositeMesh(pm, 3)
end

# @testset "2D" begin
# # given dispersion function accept a k in cartesian
# # goal is to find k_F at direction specified by angle
# dispersion(k) = dot(k, k) - 1.0

# # 2d
# N = 10
# bound = [-π, π]
# theta = SimpleGrid.Uniform(bound, N; isperiodic=true)

# k_F_previous = 0.0
# for θ in theta
# f(k) = dispersion(BZMeshes._polar2cart(Polar(k, θ)))
# k_F = find_zero(f, k_F_previous)
# @test k_F ≈ 1.0
# @test find_kFermi(dispersion, θ; kinit=k_F_previous) ≈ 1.0
# k_F_previous = k_F
# end

# # grids = kF_densed_kgrids(dispersion=dispersion, anglemesh=theta,
# # bound=[0.0, 2.0])
# # cm = CompositeMesh(theta, grids)
# DIM = 2
# lattice = Matrix([1.0 0; 0 1]')
# br = BZMeshes.Cell(lattice=lattice)

# pm = PolarMesh(dispersion=dispersion, anglemesh=theta, cell=br, kmax=2.0)
# @test AbstractMeshes.volume(pm) ≈ 4π

# end

# @testset "3D" begin
# dispersion(k) = dot(k, k) - 1.0

# N = 6
# bound = [-π, π]
# phi = SimpleGrid.Uniform(bound, N; isperiodic=true)

# N = 4
# bound = [-π / 2, π / 2]
# theta = SimpleGrid.Uniform(bound, N; isperiodic=true)

# am = CompositeMesh(phi, [theta for i in 1:length(phi)])

# k_F_previous = 0.0
# for ap in am
# f(k) = dispersion(BZMeshes._spherical2cart(Spherical(k, ap...)))
# k_F = find_zero(f, k_F_previous)
# @test k_F ≈ 1.0
# @test find_kFermi(dispersion, ap; kinit=k_F_previous) ≈ 1.0
# k_F_previous = k_F
# end

# DIM = 3
# lattice = Matrix([1.0 1.0 0; 1 0 1; 0 1 1]')
# br = BZMeshes.Cell(lattice=lattice)

# pm = PolarMesh(dispersion=dispersion, anglemesh=am, cell=br, kmax=2.0)
# @test AbstractMeshes.volume(pm) ≈ 32π / 3

# end

# @testset "Radial rescale" begin
# @testset "RescaledGrid" begin
# rmax = 2.0
# N = 10
# rgrid = SimpleG.Uniform([0.0, rmax^2], N; gpbound=[rmax^2 / 2N, rmax^2 * (1 - 1 / 2N)])
# rrg = radial_rescale(grid=rgrid, DIM=2)
# println(rrg.grid)
# bound = [-π, π]
# theta = SimpleGrid.Uniform(bound, N; gpbound=[π * (1 / N - 1), π * (1 - 1 / N)], isperiodic=true)
pm2 = BZMeshes.CompositePolarMesh(dispersion=dispersion,
anglemesh=theta, cell=br, basegridtype=:uniform,
kmax=π * sqrt(2.4), Nloggrid=10, Nbasegrid=6, minterval=0.01T, N=6)

# # for (i, p) in enumerate(theta)
# # println(volume(theta, i))
# # end

# DIM = 2
# lattice = Matrix([1.0 0; 0 1]')
# br = BZMeshes.Cell(lattice=lattice)
# bzmesh = PolarMesh(br, CompositeMesh(theta, [rrg for i in 1:length(theta)]))
# for (i, p) in enumerate(bzmesh)
# println(p, AbstractMeshes.volume(bzmesh, i))
# end
# end
println("N=$(length(pm)), ", int_mesh(pm))
println("N=$(length(pm2)), ", int_mesh(pm2))

# end
end
Nuni = 1000
um = BZMeshes.UniformBZMesh(cell=br, size=(Nuni, Nuni))
um2 = BZMeshes.UniformBZMesh(cell=br, size=(2Nuni, 2Nuni))
println("N=$(length(um)), ", int_mesh(um))
println("N=$(length(um2)), ", int_mesh(um2))
10 changes: 10 additions & 0 deletions src/BaseMesh.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,16 @@ Here we assume periodic boundary condition so for all case it's the same.
AbstractMeshes.volume(mesh::AbstractUniformMesh) = cell_volume(mesh)
AbstractMeshes.volume(mesh::AbstractUniformMesh, i) = cell_volume(mesh) / length(mesh)

AbstractMeshes.interp(data, mesh::AbstractUniformMesh, x) = data[locate(mesh, x)]

function AbstractMeshes.integrate(data, mesh::AbstractUniformMesh)
result = 0.0
for i in 1:length(mesh)
result += data[i] * volume(mesh, i)
end
return result
end

struct HasLattice <: AbstractMeshes.LatticeStyle end # has mesh.lattice and mesh.inv_lattice
AbstractMeshes.lattice_vector(::HasLattice, mesh) = mesh.lattice
AbstractMeshes.inv_lattice_vector(::HasLattice, mesh) = mesh.inv_lattice
Expand Down
4 changes: 1 addition & 3 deletions src/PolarMeshes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ function AbstractMeshes.integrate(data, mesh::PolarMesh{T,3,MT}) where {T,MT}
weighteddata = similar(data)
for pi in 1:length(mesh)
r, θ, ϕ = _extract(mesh[AngularCoords, pi])
weighteddata[pi] = r^2 * sin(θ) * data[pi]
weighteddata[pi] = r^2 * cos(θ) * data[pi]
end
return integrate(weighteddata, mesh.mesh)
end
Expand Down Expand Up @@ -270,7 +270,6 @@ function PolarMesh(; dispersion, anglemesh, cell, kmax,

DIM = size(cell.lattice, 1)
bound = [0.0, kmax]
println(typeof(dispersion))
grids = kF_densed_kgrids(; dispersion=dispersion, anglemesh=anglemesh, bound=bound, DIM=DIM, kwargs...)
cm = ProdMesh(grids, anglemesh)
pm = PolarMesh(cell, cm)
Expand All @@ -282,7 +281,6 @@ function CompositePolarMesh(; dispersion, anglemesh, cell, kmax, N,

DIM = size(cell.lattice, 1)
bound = [0.0, kmax]
println(typeof(dispersion))
grids = kF_densed_kgrids(; dispersion=dispersion, anglemesh=anglemesh, bound=bound, DIM=DIM, kwargs...)
prm = ProdMesh(grids, anglemesh)
cm = CompositeMesh(prm, N)
Expand Down
38 changes: 38 additions & 0 deletions test/BZMeshes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,6 +231,44 @@

end

@testset "3D CompositePolarMesh" begin
dispersion(k) = dot(k, k) - 1.0

N = 6
bound = [-π, π]
phi = SimpleGrid.Uniform(bound, N)

N = 4
bound = [-π / 2, π / 2]
theta = SimpleGrid.Uniform(bound, N)

am = ProdMesh([theta for i in 1:length(phi)], phi)

DIM = 3
lattice = Matrix([1.0 1.0 0; 1 0 1; 0 1 1]')
br = BZMeshes.Cell(lattice=lattice)

pm = CompositePolarMesh(dispersion=dispersion, anglemesh=am, cell=br, kmax=2.0, N=4)
@test AbstractMeshes.volume(pm) ≈ 32π / 3

data = zeros(size(pm))
for (pi, p) in enumerate(pm)
data[pi] = dispersion(p)
end

testN = 10
for i in 1:testN
r, ϕ, θ = rand(rng) * 2.0, (rand(rng) * 2 - 1) * π, (rand(rng) * 2 - 1) * π / 2
p = Spherical(r, θ, ϕ)
x = BZMeshes._cartesianize(p)
@test isapprox(AbstractMeshes.interp(data, pm, x), dispersion(x), rtol=1e-4)
@test isapprox(AbstractMeshes.interp(data, pm, p), dispersion(x), rtol=1e-4)
end

@test isapprox(AbstractMeshes.integrate(data, pm), 4π * 56 / 15, rtol=1e-4)

end

@testset "Radial rescale" begin
@testset "RescaledGrid" begin
rmax = 2.0
Expand Down

0 comments on commit 27b5227

Please # to comment.