Skip to content

Commit

Permalink
Add edge reaction term to implement joule heat (#61)
Browse files Browse the repository at this point in the history 
* implement Joule heating example
* add edge reaction test example
* test function handling for edge reactions
  • Loading branch information
@j-fu
j-fu authored Feb 24, 2023
1 parent 7446c8e commit f50f307
Show file tree
Hide file tree
Showing 2 changed files with 204 additions and 0 deletions.
178 changes: 178 additions & 0 deletions examples/Example002_EdgeReaction.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,178 @@
# # 002: check edge reaction
# ([source code](SOURCE_URL))

module Example002_EdgeReaction

using Printf
using VoronoiFVM
using ExtendableGrids
using ExtendableSparse
using GridVisualize
using LinearAlgebra
using SimplexGridFactory
using Triangulate

function main(; nref = 0, dim=2, Plotter = nothing, verbose = "and", case=:compare_max)



X=0:0.25*2.0^-nref:1
i0::Int=0
i1::Int=0
if dim==1
grid=simplexgrid(X)
i0=1
i1=2
elseif dim==2
b=SimplexGridBuilder(Generator=Triangulate)
p00=point!(b,0,0)
p10=point!(b,1,0)
p11=point!(b,1,1)
p01=point!(b,0,1)
pxx=point!(b,0.3,0.3)

facetregion!(b,1)
facet!(b,p00,p10)
facetregion!(b,2)
facet!(b,p10,p11)
facetregion!(b,3)
facet!(b,p11,p01)
facetregion!(b,4)
facet!(b,p01,p00)
grid=simplexgrid(b;maxvolume=0.01*4.0^(-nref))
i0=1
i1=3
elseif dim==3
grid=simplexgrid(X,X,X)
i0=5
i1=6
end


function storage!(y,u,node)
y[1]=u[1]
end

function flux!(y,u,edge)
y[1]= u[1,1]-u[1,2]
end

# Three ways to give a constant reaction term. As a consequence,
# these need to yield the same solution.
# 1: classical node reaction, multiplied by control volume size
function reaction!(y,u,node)
y[1]=-1
end

# 2: Edge reaction. Here we give it as a constant, and wie need
# to turn the multiplication with σ/h into a multiplication with the
# half diamond volume.
#
# Half diamond volume calculation
# /|\
# / | \
# / |s \
# -------
# h
# A=s*h/2d . Our formfactor: σ=s/h => A=σ*h^2
# - make transfer area to volume
#
# τ=1/h v= s*h/2d = σ*h^2/2d
#
function edgereaction!(y,u,edge)
h=meas(edge)
y[1]=-1*h^2/(2*dim)
end

#
# 3: "Joule heat: |∇ϕ|^2=1 after 3.17 in Bradji/Herbin
# Here we divide twice by "h" to get the constant squared gradient.
# The multiplication with dim in 3.17 compensates the division
# we had before
ϕ=grid[Coordinates][1,:]

function edgereaction2!(y,u,edge)
ϕK=ϕ[edge.node[1]]
ϕL=ϕ[edge.node[2]]
y[1]=-(ϕK-ϕL)*(ϕK-ϕL)/2
end

if case==:compare_max

function bcondition!(y,u,node)
boundary_dirichlet!(y,u,node,species=1,region=1, value=0)
boundary_dirichlet!(y,u,node,species=1,region=2, value=0)
boundary_dirichlet!(y,u,node,species=1,region=3, value=0)
boundary_dirichlet!(y,u,node,species=1,region=4, value=0)
boundary_dirichlet!(y,u,node,species=1,region=5, value=0)
boundary_dirichlet!(y,u,node,species=1,region=6, value=0)
end

sys_noderea=VoronoiFVM.System(grid; bcondition=bcondition!, flux=flux!, reaction=reaction!, storage=storage!,species=[1], is_linear=true)
sys_edgerea=VoronoiFVM.System(grid; bcondition=bcondition!, flux=flux!, edgereaction=edgereaction!, storage=storage!,species=[1], is_linear=true)
sys_edgerea2=VoronoiFVM.System(grid; bcondition=bcondition!, flux=flux!, edgereaction=edgereaction2!, storage=storage!,species=[1], is_linear=true)

sol_noderea=solve(sys_noderea;verbose)
sol_edgerea=solve(sys_edgerea;verbose)
sol_edgerea2=solve(sys_edgerea2;verbose)

vis=GridVisualizer(;Plotter,layout=(2,2))
scalarplot!(vis[1,1],grid,sol_noderea[1,:],title="node reaction",colormap=:hot)
scalarplot!(vis[2,1],grid,sol_edgerea[1,:],title="edgerea1",colormap=:hot)
scalarplot!(vis[1,2],grid,sol_edgerea2[1,:],title="edgerea2",colormap=:hot)

reveal(vis)
return maximum.([sol_noderea,sol_edgerea,sol_edgerea2])
end

if case==:compare_flux
function bcondition2!(y,u,node)
boundary_dirichlet!(y,u,node,species=1,region=i1, value=0)
end


sys2_noderea=VoronoiFVM.System(grid; bcondition=bcondition2!, flux=flux!, reaction=reaction!, storage=storage!,species=[1], is_linear=true)
sys2_edgerea=VoronoiFVM.System(grid; bcondition=bcondition2!, flux=flux!, edgereaction=edgereaction!, storage=storage!,species=[1], is_linear=true)
sys2_edgerea2=VoronoiFVM.System(grid; bcondition=bcondition2!, flux=flux!, edgereaction=edgereaction2!, storage=storage!,species=[1], is_linear=true)

sol2_noderea=solve(sys2_noderea;verbose)
sol2_edgerea=solve(sys2_edgerea;verbose)
sol2_edgerea2=solve(sys2_edgerea2;verbose)

tfac2_noderea=TestFunctionFactory(sys2_noderea)
tfc2_noderea=testfunction(tfac2_noderea,[i0],[i1])

tfac2_edgerea=TestFunctionFactory(sys2_edgerea)
tfc2_edgerea=testfunction(tfac2_edgerea,[i0],[i1])

tfac2_edgerea2=TestFunctionFactory(sys2_edgerea2)
tfc2_edgerea2=testfunction(tfac2_edgerea2,[i0],[i1])

vis=GridVisualizer(;Plotter,layout=(2,2))
scalarplot!(vis[1,1],grid,sol2_noderea[1,:],title="node reaction",colormap=:hot)
scalarplot!(vis[2,1],grid,sol2_edgerea[1,:],title="edgerea1",colormap=:hot)
scalarplot!(vis[1,2],grid,sol2_edgerea2[1,:],title="edgerea2",colormap=:hot)
reveal(vis)

I_noderea=integrate(sys2_noderea,tfc2_noderea,sol2_noderea)
I_edgerea=integrate(sys2_edgerea,tfc2_edgerea,sol2_edgerea)
I_edgerea2=integrate(sys2_edgerea2,tfc2_edgerea2,sol2_edgerea2)

return I_noderea,I_edgerea, I_edgerea2
end
end

function test()
res=fill(false,3)
for dim in 1:3
result_max=main(case=:compare_max)
result_flux=main(case=:compare_flux)
res[dim]=isapprox(result_max[1],result_max[2],atol=1.0e-6) &&
isapprox(result_max[1],result_max[3],atol=1.0e-3) &&
isapprox(result_flux[1],result_flux[2],atol=1.0e-10) &&
isapprox(result_flux[1],result_flux[3],atol=1.0e-10)
end
all(a->a,res)
end

end
26 changes: 26 additions & 0 deletions src/vfvm_testfunctions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,7 @@ function integrate(system::AbstractSystem, tf, U::AbstractMatrix{Tv},

src_eval = ResEvaluator(physics, :source, UK, node, nspecies + nparams)
rea_eval = ResEvaluator(physics, :reaction, UK, node, nspecies + nparams)
erea_eval = ResEvaluator(physics, :edgereaction, UK, edge, nspecies + nparams)
stor_eval = ResEvaluator(physics, :storage, UK, node, nspecies + nparams)
storold_eval = ResEvaluator(physics, :storage, UKold, node, nspecies + nparams)
flux_eval = ResEvaluator(physics, :flux, UKL, edge, nspecies + nparams)
Expand All @@ -140,6 +141,17 @@ function integrate(system::AbstractSystem, tf, U::AbstractMatrix{Tv},
integral[ispec] += system.celledgefactors[iedge, icell] * flux[ispec] * (tf[edge.node[1]] - tf[edge.node[2]])
end
assemble_res(edge, system, asm_res)

if isnontrivial(erea_eval)
evaluate!(erea_eval, UKL)
erea = res(erea_eval)

function easm_res(idofK, idofL, ispec)
integral[ispec] += system.celledgefactors[iedge, icell] * erea[ispec] * (tf[edge.node[1]] + tf[edge.node[2]])
end
assemble_res(edge, system, easm_res)

end
end

for inode = 1:num_nodes(geom)
Expand Down Expand Up @@ -201,6 +213,7 @@ function integrate_stdy(system::AbstractSystem, tf::Vector{Tv}, U::AbstractArray

src_eval = ResEvaluator(physics, :source, UK, node, nspecies)
rea_eval = ResEvaluator(physics, :reaction, UK, node, nspecies)
erea_eval = ResEvaluator(physics, :edgereaction, UK, node, nspecies)
flux_eval = ResEvaluator(physics, :flux, UKL, edge, nspecies)

for icell = 1:num_cells(grid)
Expand All @@ -216,6 +229,19 @@ function integrate_stdy(system::AbstractSystem, tf::Vector{Tv}, U::AbstractArray
integral[ispec] += system.celledgefactors[iedge, icell] * flux[ispec] * (tf[edge.node[1]] - tf[edge.node[2]])
end
assemble_res(edge, system, asm_res)

if isnontrivial(erea_eval)
evaluate!(erea_eval, UKL)
erea = res(erea_eval)

function easm_res(idofK, idofL, ispec)
integral[ispec] += system.celledgefactors[iedge, icell] * erea[ispec] * (tf[edge.node[1]] + tf[edge.node[2]])
end
assemble_res(edge, system, easm_res)

end


end

for inode = 1:num_nodes(geom)
Expand Down

2 comments on commit f50f307

@j-fu
Copy link
Member Author

@j-fu j-fu commented on f50f307 Feb 24, 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/78432

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.1.0 -m "<description of version>" f50f3076b456d4845e6a691e43fa9d9eb354cbda
git push origin v1.1.0

Please sign in to comment.