Skip to content

Commit

Permalink
_complete system in postprocessing methods
Browse files Browse the repository at this point in the history
  • Loading branch information
@j-fu
j-fu committed Dec 11, 2023
1 parent 9b1880a commit 8daba8c
Show file tree
Hide file tree
Showing 4 changed files with 137 additions and 141 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "VoronoiFVM"
uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
authors = ["Jürgen Fuhrmann <juergen.fuhrmann@wias-berlin.de>", "Dilara Abdel", "Jan Weidner", "Alexander Seiler", "Patricio Farrell", "Matthias Liero"]
version = "1.15.0"
version = "1.15.1"

[deps]
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
Expand Down
141 changes: 71 additions & 70 deletions src/vfvm_postprocess.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ The result is an `nspec x nregion` vector.
"""
function integrate(system::AbstractSystem{Tv, Tc, Ti, Tm}, F::Function, U::AbstractMatrix{Tu};
boundary = false) where {Tu, Tv, Tc, Ti, Tm}
_complete!(system)
grid = system.grid
data = system.physics.data
nspecies = num_species(system)
Expand All @@ -35,8 +36,8 @@ function integrate(system::AbstractSystem{Tv, Tc, Ti, Tm}, F::Function, U::Abstr
integral = zeros(Tu, nspecies, nbfaceregions)

for item in nodebatch(system.boundary_assembly_data)
for ibnode in noderange(system.boundary_assembly_data,item)
_fill!(bnode,system.boundary_assembly_data,ibnode,item)
for ibnode in noderange(system.boundary_assembly_data, item)
_fill!(bnode, system.boundary_assembly_data, ibnode, item)
res .= zero(Tv)
@views F(rhs(bnode, res), unknowns(bnode, U[:, bnode.index]), bnodeparams...)
asm_res(idof, ispec) = integral[ispec, bnode.region] += bnode.fac * res[ispec]
Expand All @@ -55,8 +56,8 @@ function integrate(system::AbstractSystem{Tv, Tc, Ti, Tm}, F::Function, U::Abstr
ncellregions = maximum(cellregions)
integral = zeros(Tu, nspecies, ncellregions)
for item in nodebatch(system.assembly_data)
for inode in noderange(system.assembly_data,item)
_fill!(node,system.assembly_data,inode,item)
for inode in noderange(system.assembly_data, item)
_fill!(node, system.assembly_data, inode, item)
res .= zero(Tv)
@views F(rhs(node, res), unknowns(node, U[:, node.index]), nodeparams...)
asm_res(idof, ispec) = integral[ispec, node.region] += node.fac * res[ispec]
Expand All @@ -77,8 +78,9 @@ The result is an `nspec x nregion` vector.
"""
function edgeintegrate(system::AbstractSystem{Tv, Tc, Ti, Tm}, F::Function, U::AbstractMatrix{Tu};
boundary = false) where {Tu, Tv, Tc, Ti, Tm}
_complete!(system)
grid = system.grid
dim=dim_space(grid)
dim = dim_space(grid)
data = system.physics.data
nspecies = num_species(system)
res = zeros(Tu, nspecies)
Expand All @@ -98,16 +100,16 @@ function edgeintegrate(system::AbstractSystem{Tv, Tc, Ti, Tm}, F::Function, U::A
ncellregions = maximum(cellregions)
integral = zeros(Tu, nspecies, ncellregions)
for item in edgebatch(system.assembly_data)
for iedge in edgerange(system.assembly_data,item)
_fill!(edge,system.assembly_data,iedge,item)
for iedge in edgerange(system.assembly_data, item)
_fill!(edge, system.assembly_data, iedge, item)
@views UKL[1:nspecies] .= U[:, edge.node[1]]
@views UKL[(nspecies+1):(2*nspecies)] .= U[:, edge.node[2]]
@views UKL[(nspecies + 1):(2 * nspecies)] .= U[:, edge.node[2]]
res .= zero(Tv)
@views F(rhs(edge, res), unknowns(edge,UKL), edgeparams...)
@views F(rhs(edge, res), unknowns(edge, UKL), edgeparams...)
function asm_res(idofK, idofL, ispec)
h=meas(edge)
h = meas(edge)
# This corresponds to the multiplication with the diamond volume.
integral[ispec,edge.region] += h^2*edge.fac*res[ispec]/dim
integral[ispec, edge.region] += h^2 * edge.fac * res[ispec] / dim
end
assemble_res(edge, system, asm_res)
end
Expand All @@ -116,7 +118,6 @@ function edgeintegrate(system::AbstractSystem{Tv, Tc, Ti, Tm}, F::Function, U::A
return SolutionIntegral(integral)
end


############################################################################
"""
$(SIGNATURES)
Expand Down Expand Up @@ -147,6 +148,7 @@ corners in the code. Please see any potential boundary artifacts as a manifestat
of this issue and report them.
"""
function nodeflux(system::AbstractSystem{Tv, Tc, Ti, Tm}, U::AbstractArray{Tu, 2}) where {Tu, Tv, Tc, Ti, Tm}
_complete!(system)
grid = system.grid
dim = dim_space(grid)
nnodes = num_nodes(grid)
Expand All @@ -161,14 +163,13 @@ function nodeflux(system::AbstractSystem{Tv, Tc, Ti, Tm}, U::AbstractArray{Tu, 2
edge = Edge(system)
node = Node(system)


# !!! TODO Parameter handling here
UKL = Array{Tu, 1}(undef, 2 * nspecies)
flux_eval = ResEvaluator(system.physics, :flux, UKL, edge, nspecies)

for item in edgebatch(system.assembly_data)
for iedge in edgerange(system.assembly_data,item)
_fill!(edge,system.assembly_data,iedge,item)
for iedge in edgerange(system.assembly_data, item)
_fill!(edge, system.assembly_data, iedge, item)
K = edge.node[1]
L = edge.node[2]
@views UKL[1:nspecies] .= U[:, edge.node[1]]
Expand All @@ -182,171 +183,171 @@ function nodeflux(system::AbstractSystem{Tv, Tc, Ti, Tm}, U::AbstractArray{Tu, 2
assemble_res(edge, system, asm_res)
end
end

for item in nodebatch(system.assembly_data)
for inode in noderange(system.assembly_data,item)
_fill!(node,system.assembly_data,inode,item)
for inode in noderange(system.assembly_data, item)
_fill!(node, system.assembly_data, inode, item)
nodevol[node.index] += node.fac
end
end

for inode = 1:nnodes
@views nodeflux[:, :, inode] /= nodevol[inode]
end
nodeflux
end


#####################################################################################
"""
$(SIGNATURES)
Calculate Euklidean norm of the degree of freedom vector.
"""
function LinearAlgebra.norm(system::AbstractSystem, u, p::Number=2) end
function LinearAlgebra.norm(system::AbstractSystem, u, p::Number = 2) end

function LinearAlgebra.norm(system::DenseSystem, u, p::Number=2)
function LinearAlgebra.norm(system::DenseSystem, u, p::Number = 2)
_initialize_inactive_dof!(u, system)
_complete!(system)
norm(u, p)
end


LinearAlgebra.norm(system::SparseSystem, u::SparseSolutionArray, p::Number=2) = LinearAlgebra.norm(u.node_dof.nzval, p)
LinearAlgebra.norm(system::SparseSystem, u::SparseSolutionArray, p::Number = 2) = LinearAlgebra.norm(u.node_dof.nzval, p)

"""
$(SIGNATURES)
Calculate weighted discrete ``L^p`` norm of a solution vector.
"""
function lpnorm(sys,u,p,species_weights=ones(num_species(sys)))
nspec=num_species(sys)
II=integrate(sys, (y,u,node,data=nothing)->y.=u.^p,u)
(sum([sum(II[i,:]) for i=1:nspec].*species_weights))^(1/p)
function lpnorm(sys, u, p, species_weights = ones(num_species(sys)))
_complete!(sys)
nspec = num_species(sys)
II = integrate(sys, (y, u, node, data = nothing) -> y .= u .^ p, u)
(sum([sum(II[i, :]) for i = 1:nspec] .* species_weights))^(1 / p)
end

"""
$(SIGNATURES)
Calculate weigthed discrete ``L^2(\\Omega)`` norm of a solution vector.
"""
function l2norm(sys,u,species_weights=ones(num_species(sys)))
lpnorm(sys,u,2,species_weights)
function l2norm(sys, u, species_weights = ones(num_species(sys)))
_complete!(sys)
lpnorm(sys, u, 2, species_weights)
end

"""
$(SIGNATURES)
Calculate weighted discrete ``W^{1,p}(\\Omega)`` seminorm of a solution vector.
"""
function w1pseminorm(sys,u,p,species_weights=ones(num_species(sys)))
nspec=num_species(sys)
dim=dim_space(sys.grid)
function f(y,u,edge,data=nothing)
h=meas(edge)
for ispec=1:nspec
y[ispec]=dim*((u[ispec,1]-u[ispec,2])/h)^p
function w1pseminorm(sys, u, p, species_weights = ones(num_species(sys)))
_complete!(sys)
nspec = num_species(sys)
dim = dim_space(sys.grid)
function f(y, u, edge, data = nothing)
h = meas(edge)
for ispec = 1:nspec
y[ispec] = dim * ((u[ispec, 1] - u[ispec, 2]) / h)^p
end
end
II=edgeintegrate(sys,f,u)
(sum([sum(II[i,:]) for i=1:nspec].*species_weights))^(1/p)
II = edgeintegrate(sys, f, u)
(sum([sum(II[i, :]) for i = 1:nspec] .* species_weights))^(1 / p)
end

"""
$(SIGNATURES)
Calculate weighted discrete ``H^1(\\Omega)`` seminorm of a solution vector.
"""
function h1seminorm(sys,u,species_weights=ones(num_species(sys)))
w1pseminorm(sys,u,2,species_weights)
function h1seminorm(sys, u, species_weights = ones(num_species(sys)))
_complete!(sys)
w1pseminorm(sys, u, 2, species_weights)
end

"""
$(SIGNATURES)
Calculate weighted discrete ``W^{1,p}(\\Omega)`` norm of a solution vector.
"""
function w1pnorm(sys,u,p,species_weights=ones(num_species(sys)))
lpnorm(sys,u,p,species_weights)+w1pseminorm(sys,u,p,species_weights)
function w1pnorm(sys, u, p, species_weights = ones(num_species(sys)))
_complete!(sys)
lpnorm(sys, u, p, species_weights) + w1pseminorm(sys, u, p, species_weights)
end

"""
$(SIGNATURES)
Calculate weighted discrete ``H^1(\\Omega)`` norm of a solution vector.
"""
function h1norm(sys,u,species_weights=ones(num_species(sys)))
w1pnorm(sys,u,2,species_weights)
function h1norm(sys, u, species_weights = ones(num_species(sys)))
_complete!(sys)
w1pnorm(sys, u, 2, species_weights)
end

function _bochnernorm(sys,u::TransientSolution,p,species_weights,spatialnorm::F) where F
n=length(u.t)
nrm=spatialnorm(sys,u.u[1],p,species_weights)^p*(u.t[2]-u.t[1])/2
for i=2:n-1
nrm+=spatialnorm(sys,u.u[i],p,species_weights)^p*(u.t[i+1]-u.t[i-1])/2
function _bochnernorm(sys, u::TransientSolution, p, species_weights, spatialnorm::F) where {F}
_complete!(sys)
n = length(u.t)
nrm = spatialnorm(sys, u.u[1], p, species_weights)^p * (u.t[2] - u.t[1]) / 2
for i = 2:(n - 1)
nrm += spatialnorm(sys, u.u[i], p, species_weights)^p * (u.t[i + 1] - u.t[i - 1]) / 2
end
nrm+=spatialnorm(sys,u.u[n],p,species_weights)^p*(u.t[end]-u.t[end-1])/2
nrm^(1/p)
nrm += spatialnorm(sys, u.u[n], p, species_weights)^p * (u.t[end] - u.t[end - 1]) / 2
nrm^(1 / p)
end

"""
$(SIGNATURES)
Calculate weighted discrete ``L^p([0,T];W^{1,p}(\\Omega))`` norm of a transient solution.
"""
function lpw1pnorm(sys,u::TransientSolution,p,species_weights=ones(num_species(sys)))
_bochnernorm(sys,u,p,species_weights,w1pnorm)
function lpw1pnorm(sys, u::TransientSolution, p, species_weights = ones(num_species(sys)))
_bochnernorm(sys, u, p, species_weights, w1pnorm)
end


"""
$(SIGNATURES)
Calculate weighted discrete ``L^p([0,T];W^{1,p}(\\Omega))`` seminorm of a transient solution.
"""
function lpw1pseminorm(sys,u::TransientSolution,p,species_weights=ones(num_species(sys)))
_bochnernorm(sys,u,p,species_weights,w1pseminorm)
function lpw1pseminorm(sys, u::TransientSolution, p, species_weights = ones(num_species(sys)))
_bochnernorm(sys, u, p, species_weights, w1pseminorm)
end



"""
$(SIGNATURES)
Calculate weighted discrete ``L^2([0,T];H^1(\\Omega))`` seminorm of a transient solution.
"""
function l2h1seminorm(sys,u::TransientSolution,species_weights=ones(num_species(sys)))
lpw1pseminorm(sys,u,2,species_weights)
function l2h1seminorm(sys, u::TransientSolution, species_weights = ones(num_species(sys)))
lpw1pseminorm(sys, u, 2, species_weights)
end


"""
$(SIGNATURES)
Calculate weighted discrete ``L^2([0,T];H^1(\\Omega))`` norm of a transient solution.
"""
function l2h1norm(sys,u::TransientSolution,species_weights=ones(num_species(sys)))
lpw1pnorm(sys,u,2,species_weights)
function l2h1norm(sys, u::TransientSolution, species_weights = ones(num_species(sys)))
lpw1pnorm(sys, u, 2, species_weights)
end


"""
$(SIGNATURES)
Calculate volumes of Voronoi cells.
"""
function nodevolumes(system)
_complete!(system)
if isnothing(system.assembly_data)
update_grid!(system)
end
node= Node(system)
node = Node(system)
nodevol = zeros(num_nodes(system.grid))
for item in nodebatch(system.assembly_data)
for inode in noderange(system.assembly_data,item)
_fill!(node,system.assembly_data,inode,item)
for inode in noderange(system.assembly_data, item)
_fill!(node, system.assembly_data, inode, item)
nodevol[node.index] += node.fac
end
end
nodevol
end


17 changes: 13 additions & 4 deletions src/vfvm_solver.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -352,6 +352,7 @@ function CommonSolve.solve(inival,
while !solved
solved = true
forced = false
errored = false
try # check for non-converging newton
λ = λ0 + Δλ
control.pre(solution, λ)
Expand All @@ -377,15 +378,14 @@ function CommonSolve.solve(inival,
params = params,
kwargs...,)
catch err
if transient
err = "Problem at $(λstr)=$(λ|>rd), Δ$(λstr)=$(Δλ|>rd):\n$(err)"
end
err = "Problem at $(λstr)=$(λ|>rd), Δ$(λstr)=$(Δλ|>rd):\n$(err)"
if (control.handle_exceptions)
_print_error(err, stacktrace(catch_backtrace()))
else
rethrow(err)
end
solved = false
errored = true
end
if solved
Δu = control.delta(system, solution, oldsolution, λ, Δλ)
Expand All @@ -406,12 +406,21 @@ function CommonSolve.solve(inival,
throw(ErrorException(err))
end
break # give up lowering stepsize, break out if "while !solved" loop
else
elseif !errored
if doprint(control, 'e')
println("[e]volution: forced first timestep: Δu/Δu_opt=$(Δu/Δu_opt|>rd)")
end
forced = true
solved = true
else
err = "Convergence problem in first timestep"
if control.handle_exceptions
@warn err
else
throw(ErrorException(err))
end

break
end
else
# reduce time step
Expand Down
Loading

2 comments on commit 8daba8c

@j-fu
Copy link
Owner Author

@j-fu j-fu commented on 8daba8c Dec 11, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/96870

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v1.15.1 -m "<description of version>" 8daba8cd1bcefeddcc43bd98673a1bb141aafb67
git push origin v1.15.1

Please sign in to comment.