|
| 1 | +using Lehmann |
| 2 | +using StaticArrays, Printf |
| 3 | +using CompositeGrids |
| 4 | +using LinearAlgebra |
| 5 | +using DelimitedFiles |
| 6 | +#const Float = BigFloat #FQR.Float |
| 7 | +#const Double =BigFloat #FQR.Double |
| 8 | +#const DotF = BigFloat |
| 9 | +#const Tiny = DotF(1e-5) |
| 10 | + |
| 11 | +struct TauGrid{Float} <: FQR.Grid |
| 12 | + tau::Float # actual location of the grid point |
| 13 | + coord::Int # integer coordinate of the grid point on the fine meshes |
| 14 | + vec::Vector{Float} |
| 15 | +end |
| 16 | + |
| 17 | +Base.show(io::IO, grid::TauGrid) = print(io, "τ = ($(@sprintf("%12.4f", grid.tau[1])))") |
| 18 | + |
| 19 | +struct TauFineMesh{Float} <: FQR.FineMesh |
| 20 | + symmetry::Int # symmetrize (omega1, omega2) <-> (omega2, omega1) |
| 21 | + candidates::Vector{TauGrid{Float}} # vector of grid points |
| 22 | + selected::Vector{Bool} |
| 23 | + residual::Vector{Float} |
| 24 | + |
| 25 | + ## for frequency mesh only ### |
| 26 | + fineGrid::Vector{Float} # fine grid for each dimension |
| 27 | + |
| 28 | + function TauFineMesh(U::AbstractMatrix{Float}, finegrid; sym=1) where {Float} |
| 29 | + # initialize the residual on fineGrid with <g, g> |
| 30 | + #_finegrid = Float.(τChebyGrid(Λ)) |
| 31 | + # separationTest(_finegrid) |
| 32 | + mesh = new{Float}(sym, [], [], [], finegrid) |
| 33 | + |
| 34 | + for (xi, x) in enumerate(finegrid) |
| 35 | + coord = xi |
| 36 | + #if irreducible(coord, sym, Nfine) # if grid point is in the reducible zone, then skip residual initalization |
| 37 | + vec = U[xi,:] |
| 38 | + g = TauGrid(x, coord, vec) |
| 39 | + push!(mesh.candidates, g) |
| 40 | + push!(mesh.residual, FQR.dot(mesh, g, g)) |
| 41 | + push!(mesh.selected, false) |
| 42 | + #end |
| 43 | + end |
| 44 | + |
| 45 | + println("fine mesh initialized.") |
| 46 | + return mesh |
| 47 | + end |
| 48 | +end |
| 49 | + |
| 50 | + |
| 51 | + |
| 52 | +# function τChebyGrid(Λ, degree=24, print = true) |
| 53 | +# npt = Int(ceil(log(Λ) / log(2.0))) - 2 # subintervals on [0,1/2] in tau space (# subintervals on [0,1] is 2*npt) |
| 54 | + |
| 55 | +# pbpt = zeros(Float, 2npt + 1) |
| 56 | +# pbpt[1] = 0.0 |
| 57 | +# for i = 1:npt |
| 58 | +# pbpt[i+1] = 1.0 / 2^(npt - i + 1) |
| 59 | +# end |
| 60 | +# pbpt[npt+2:2npt+1] = 1 .- pbpt[npt:-1:1] |
| 61 | +# #println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[end])]") |
| 62 | +# finegrid = Lehmann.Discrete.CompositeChebyshevGrid(degree, pbpt).grid |
| 63 | +# #println("$(finegrid)\n")#[1:(length(finegrid)÷2+1)]) |
| 64 | +# #println("$(finegrid+reverse(finegrid))\n") |
| 65 | +# return finegrid |
| 66 | +# end |
| 67 | + |
| 68 | +""" |
| 69 | +composite expoential grid |
| 70 | +""" |
| 71 | +function fine_τGrid(Λ::Float,degree,ratio::Float) where {Float} |
| 72 | + ############## use composite grid ############################################# |
| 73 | + # Generating a log densed composite grid with LogDensedGrid() |
| 74 | + npo = Int(ceil(log(Λ) / log(ratio))) - 2 # subintervals on [0,1/2] in tau space (# subintervals on [0,1] is 2*npt) |
| 75 | + grid = CompositeGrid.LogDensedGrid( |
| 76 | + :cheb,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb |
| 77 | + [0.0, 1.0],# The grid is defined on [0.0, β] |
| 78 | + [0.0, 1.0],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter. |
| 79 | + npo,# N of log grid |
| 80 | + 0.5 / ratio^(npo-1), # minimum interval length of log grid |
| 81 | + degree, # N of bottom layer |
| 82 | + Float |
| 83 | + ) |
| 84 | + #print(grid[1:length(grid)÷2+1]) |
| 85 | + #print(grid+reverse(grid)) |
| 86 | + # println("Composite expoential grid size: $(length(grid))") |
| 87 | + println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[end])]") |
| 88 | + return grid |
| 89 | + |
| 90 | + ############# DLR based fine grid ########################################## |
| 91 | + # dlr = DLRGrid(Euv=Float64(Λ), beta=1.0, rtol=Float64(rtol) / 100, isFermi=true, symmetry=:ph, rebuild=true) |
| 92 | + # # println("fine basis number: $(dlr.size)\n", dlr.ω) |
| 93 | + # degree = 4 |
| 94 | + # grid = Vector{Double}(undef, 0) |
| 95 | + # panel = Double.(dlr.τ) |
| 96 | + # for i in 1:length(panel)-1 |
| 97 | + # uniform = [panel[i] + (panel[i+1] - panel[i]) / degree * j for j in 0:degree-1] |
| 98 | + # append!(grid, uniform) |
| 99 | + # end |
| 100 | + |
| 101 | + # println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[2])]") |
| 102 | + # return grid |
| 103 | +end |
| 104 | + |
| 105 | +# """ |
| 106 | +# Test the finegrids do not overlap |
| 107 | +# """ |
| 108 | +# function separationTest(finegrid) |
| 109 | +# epsilon = eps(DotF(1)) * 10 |
| 110 | +# for (i, f) in enumerate(finegrid) |
| 111 | +# # either zero, or sufficiently large |
| 112 | +# @assert abs(f) < epsilon || abs(f) > Tiny "$i: $f should either smaller than $epsilon or larger than $Tiny" |
| 113 | +# for (j, g) in enumerate(finegrid) |
| 114 | +# # two frequencies are either the same, or well separated |
| 115 | +# @assert abs(f - g) < epsilon || abs(f - g) > Tiny "$i: $f and $j: $g should either closer than $epsilon or further than $Tiny" |
| 116 | +# fg = f + g |
| 117 | +# for (k, l) in enumerate(finegrid) |
| 118 | +# @assert abs(l - fg) < epsilon || abs(l - fg) > Tiny "$i: $f + $j: $g = $fg and $k: $l should either closer than $epsilon or further than $Tiny" |
| 119 | +# end |
| 120 | +# end |
| 121 | +# end |
| 122 | +# return |
| 123 | +# end |
| 124 | + |
| 125 | +# function irreducible(coord, symmetry,length) |
| 126 | +# @assert iseven(length) "The fineGrid should have even number of points" |
| 127 | +# if symmetry == 0 |
| 128 | +# return true |
| 129 | +# else |
| 130 | +# return coord<length÷2+1 |
| 131 | +# end |
| 132 | +# end |
| 133 | + |
| 134 | +# function FQR.irreducible(grid::TauGrid) |
| 135 | +# return irreducible(grid.coord, mesh.symmetry, length(mesh.fineGrid)) |
| 136 | +# end |
| 137 | + |
| 138 | +function FQR.mirror(mesh::TauFineMesh{Float}, idx) where {Float} |
| 139 | + meshsize = length(mesh.candidates) |
| 140 | + if mesh.symmetry == 0 |
| 141 | + return [],[-1] |
| 142 | + else |
| 143 | + newgrids = TauGrid{Float}[] |
| 144 | + idxmirror = [] |
| 145 | + #coords = unique([(idx), (meshsize - idx)]) |
| 146 | + g = deepcopy(mesh.candidates[meshsize - idx+1]) |
| 147 | + #print("\n$(mesh.candidates[meshsize - idx+1].tau+mesh.candidates[idx].tau)\n") |
| 148 | + push!(newgrids, g) |
| 149 | + push!(idxmirror,meshsize - idx+1) |
| 150 | + return newgrids,idxmirror |
| 151 | + end |
| 152 | + # end |
| 153 | +end |
| 154 | + |
| 155 | + |
| 156 | +""" |
| 157 | +basis dot |
| 158 | +""" |
| 159 | +function FQR.dot(mesh, g1::TauGrid, g2::TauGrid) |
| 160 | + # println("dot: ", g1, ", ", g2) |
| 161 | + return dot(g1.vec, g2.vec) |
| 162 | +end |
| 163 | + |
| 164 | + |
| 165 | + |
| 166 | +if abspath(PROGRAM_FILE) == @__FILE__ |
| 167 | + |
| 168 | + |
| 169 | + lambda, β, rtol = 100000, 1.0,1e-8 |
| 170 | + dlr = DLRGrid(Euv=Float64(lambda), beta=β, rtol=Float64(rtol) / 100, isFermi=true, symmetry=:sym, rebuild=false) |
| 171 | + dlrfile = "basis.dat" |
| 172 | + data = readdlm(dlrfile,'\n') |
| 173 | + FreqGrid = BigFloat.(data[:,1]) |
| 174 | + mesh = TauFineMesh{BigFloat}(lambda, FreqGrid, sym=1) |
| 175 | + |
| 176 | + # KK = zeros(3, 3) |
| 177 | + # n = (2, 2) |
| 178 | + # o = (mesh.fineGrid[n[1]], mesh.fineGrid[n[2]]) |
| 179 | + # for i in 1:3 |
| 180 | + # g1 = FreqGrid{2}(i, o, n) |
| 181 | + # for j in 1:3 |
| 182 | + # g2 = FreqGrid{2}(j, o, n) |
| 183 | + # println(g1, ", ", g2) |
| 184 | + # KK[i, j] = FQR.dot(mesh, g1, g2) |
| 185 | + # end |
| 186 | + # end |
| 187 | + # display(KK) |
| 188 | + # println() |
| 189 | + |
| 190 | + basis = FQR.Basis{TauGrid, BigFloat, BigFloat}(lambda, rtol, mesh) |
| 191 | + FQR.qr!(basis, verbose=1) |
| 192 | + |
| 193 | + # lambda, rtol = 1000, 1e-8 |
| 194 | + # mesh = TauFineMesh{D}(lambda, rtol, sym=0) |
| 195 | + # basis = FQR.Basis{D,TauGrid{D}}(lambda, rtol, mesh) |
| 196 | + # @time FQR.qr!(basis, verbose=1) |
| 197 | + |
| 198 | + FQR.test(basis) |
| 199 | + |
| 200 | + mesh = basis.mesh |
| 201 | + grids = basis.grid |
| 202 | + tau_grid = [] |
| 203 | + for (i, grid) in enumerate(grids) |
| 204 | + push!(tau_grid, grid.tau) |
| 205 | + end |
| 206 | + tau_grid = sort(BigFloat.(tau_grid)) |
| 207 | + #print(tau_grid) |
| 208 | + open("basis_τ.dat", "w") do io |
| 209 | + for i in 1:length(tau_grid) |
| 210 | + println(io, tau_grid[i]) |
| 211 | + end |
| 212 | + end |
| 213 | + # Nfine = length(mesh.fineGrid) |
| 214 | + # open("finegrid.dat", "w") do io |
| 215 | + # for i in 1:Nfine |
| 216 | + # println(io, basis.mesh.fineGrid[i]) |
| 217 | + # end |
| 218 | + # end |
| 219 | + # open("residual.dat", "w") do io |
| 220 | + # # println(mesh.symmetry) |
| 221 | + # residual = zeros(Double, Nfine, Nfine) |
| 222 | + # for i in 1:length(mesh.candidates) |
| 223 | + # if mesh.candidates[i].sector == 1 |
| 224 | + # x, y = mesh.candidates[i].coord |
| 225 | + # residual[x, y] = mesh.residual[i] |
| 226 | + # # println(x, ", ", y, " -> ", length(mirror(mesh, i))) |
| 227 | + |
| 228 | + # for grid in FQR.mirror(mesh, i) |
| 229 | + # if grid.sector == 1 |
| 230 | + # xp, yp = grid.coord |
| 231 | + # residual[xp, yp] = residual[x, y] |
| 232 | + # # println(xp, ", ", yp) |
| 233 | + # end |
| 234 | + # end |
| 235 | + # end |
| 236 | + # end |
| 237 | + |
| 238 | + # for i in 1:Nfine |
| 239 | + # for j in 1:Nfine |
| 240 | + # println(io, residual[i, j]) |
| 241 | + # end |
| 242 | + # end |
| 243 | + # end |
| 244 | +end |
0 commit comments