j-fu/voronoifvm.jl

Solution of nonlinear multiphysics partial differential equation systems using the Voronoi finite volume method

voronoifvm.jl's Issues

Allow time dependency of constitutive functions

As requested in #18

Jump Conditions

Hello Jürgen Fuhrmann,

I'm very interested in your package for the simulation of the coupled transport processes in a flow battery cell. For this I'd like to use internal jump conditions: For simplicity let's consider a one-dimensional Laplace equation with Dirichlet conditions at the domain boundaries. What is the recommended way to enforce a jump in the function value (for some prescribed value) at some location x* in the domain, whereas the left and right limits of the derivative at x* are identical?

Looking through the examples it seems that this could be achieved using a DiscontinuousQuantity. However, currently it is unclear to me how to interpret the resulting discretization and achieve a prescribed jump in the solution.

Thanks for developing this very nice package.

Discontinous quantities: unexpected behavior

Setting react =0.0 in

VoronoiFVM.jl/examples/Example420_DiscontinuousQuantities.jl

Line 104 in 8a7776f

react=(u[dspec,1]-u[dspec,2])/d1

turns the dspec the DiscontinuousQuantity in continuous one:

I'd expect to see some sort of a step function -- depending on the initial condition.
@j-fu @dilaraabdel is this behavior as intended?

Drop support of Julia 1.6

I plan to drop support of Julia 1.6 for VoronoiFVM.jl and all related packages (e.g. ExtendableGrids etc.). The minimal Julia version to be supported will be 1.9 with its Extension mechanism.

Many other packages did this step, e.g. meanwhile, LinearSolve.jl (a dependency of VoronoiFVM) requires 1.10 as minimal version, which means that Pkg will resolve to older versions of LinearSolve.jl when a lower version is needed.

Probably, Julia 1.10 will become the new LTS (longtime support) release and replace 1.6 in this role.

In general, it is advisable to work with the current stable releases of Julia.

Solve transport problem with open/closed boundary

I would like to compute the steady-state of a particle density ρ, that is moved under a velocity field v and is produced at a rate s.

-div(vρ) + s = 0

The boundary condition is that particles are free to leave the simulation volume at boundaries, but no particles can enter the simulation that way. (Does this boundary condition have a name?)

I tried

using VoronoiFVM

function flux!(f,u,edge)
    v = 1.0
    h = meas(edge)
    f[1] = if v > 0
        v*u[1, 1] * h
    else
        v*u[1, 2] * h
    end
end

function source!(out, node)
    out[1] = 1
end

physics=VoronoiFVM.Physics(num_species=1, 
    flux=flux!,
    source=source!,
)

X = range(0, stop=1, length=5)
grid=VoronoiFVM.Grid(collect(X))
sys=VoronoiFVM.System(grid,physics)

ispec = 1
enable_species!(sys,ispec,[1])
boundary_neumann!(sys, ispec, 1, 0.0) # this should ensure no particles enter from the left

inival=unknowns(sys,inival=0)
solution=unknowns(sys)
solve!(solution,inival,sys)

But it errors:

ERROR: LoadError: LinearAlgebra.SingularException(0)

How to solve this problem with VoronoiFVM.jl?

Background

The above is a simplified example of an issue I face. Really I would like to simulate an ionization chamber. The equations are:

E = -grad(φ)
ε*div(E) = - ρₑ - ρ₋ + ρ₊

∂ₜρₑ =  div(μₑ(E)*E*ρₑ) + Rₑ(E,ρₑ,ρ₋,ρ₊) + sₑ
∂ₜρ₋ =  div(μ₋(E)*E*ρ₋) + R₋(E,ρₑ,ρ₋,ρ₊) + s₋
∂ₜρ₊ = -div(μ₊(E)*E*ρ₊) + R₊(E,ρₑ,ρ₋,ρ₊) + s₊

• φ is the potential
• E is the electric field
• ρₑ, ρ₋, ρ₊ are the densities of electrons, negative and positive ions
• ε is the permittivity.
• μᵢ are the mobilities of the particle species. In the case of electrons there is a strong dependence on E.
• Rᵢ are the reaction rates. For instance, an electron could attach to and O₂ molecule and become a negative ion.
  Reactions involving electrons have a dependence on E.
• sᵢ are the source terms. They capture the liberation of charge carriers due to radiation.

The unknowns are (φ, E, ρₑ, ρ₋, ρ₊ ). The other functions/constants are known. I am interested in both time evolution and steady-state (∂ₜρᵢ = 0).

Boundary conditions

The boundary decomposes into parts that are conductive Γ₁ and parts that are not Γ₂:

∂Ω = Γ₁ ∪ Γ₂

• At the conductive parts Γ₁ we have a Dirichlet condition on the potential φ
• At the non-conductive parts Γ₂ we probably want zero von Neumann condition on E

The boundary condition for the charge carriers ρᵢ is:

• At the conductive parts Γ₁ charge carriers can leave the chamber if their sign is appropriate. They can never enter the chamber.
• At the non-conductive parts Γ₂ charge carriers can neither leave nor enter.

Can I use del.f(u) (advection) term here?

Hello.
I'm studying to use this package and thank you for making the FVM model.
I'm quite a newbie for Julia so inferring model from examples and Youtube... is this right form??

Then, how can I deal with a del.f(u) term? (advection term)

Compatibility with multithreading

I wanted to look at a range of parameters on the same system, so I tried using Threads.@threads on the for loop, and I got the following error:

    nested task error: Allocations in assembly loop: cells: 208, bfaces: 0
    See the documentation of `check_allocs!` for more information

It runs just fine if I run it in series, but I thought I would try to speed up the process any way I could. Is this something VoronoiFVM cannot handle (yet), or did I fumble somewhere? Thank you!

Bug(s) in src/vfvm_assemblydata.jl

In two functions in src/vfvm_assemblydata.jl the iparam-for-loop to be moved into some other loops such that all used variables are defined.

• line 337: @inline function assemble_res_jac(edge... the iparam-for-loop has to be inside the isnodespecies(system, ispec, K) if-statement, otherwise idofL is not defined.

• line 400: @inline function assemble_res_jac(bedge... the iparam-for-loop has to be inside idofK for-loop and the isnodespecies(system, ispec, K) if-statement, otherwise idofK and idofL are not defined.

In the first function:

Precompilation

There have been issues with precompilation MacOS Intel getting stuck.
After some tests which showed no significant influence on loading time (1.10 beta, Linux) I switched off precompilation until I better understand what could be the issues and how much precompilation can help at all. See #80

hypoellictic operator

Hi,

I would like to use your package to solve an equation of the form

d/dt P = d2/dx2 P + d/dx( F(x,y) P) + d/dy (G(x,y) P)

with Dirichlet BC. However, I am struggling writing the flux function mainly because the "Laplacian" only applies to the x variable. Can you give me a hint please?

Thank you

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Try out DifferentiationInterface for sparse autodiff?

Hi @j-fu!
As of this morning's release, DifferentiationInterface is starting to look like a solid replacement for SparseDiffTools in the realm of sparse Jacobians and Hessians. Would you be interested in trying it out?
You can check out the tutorial to get a feel for the interface.

Code comments in examples

The examples need some reformatting, such that the comments inside code (e.g. inside of modules) are not becoming text (cf. fredrikekre/Literate.jl#37).
This hurts copy-pasteability a bit.

Provide function which returns solution on boundary region

For some post processing purposes, it could be beneficial to provide a function which returns the nodal values $u_k$ at designated boundary face regions.

My suggestion would be to provide something like a bvals function in vfvm_system.jl.

A quick and dirty draft could look like this:

function bvalinds(system::VoronoiFVM.AbstractSystem{Tv,Ti}, regions::AbstractVector) where {Tv,Ti}
   TA=Array{Ti}
   bvals=Array{TA}(undef,length(regions))
   grid=system.grid

   for (i,ireg) in enumerate(regions)
      catBFaces=findall(grid[BFaceRegions].==ireg)
      catBFaceLeftNodes=grid[BFaceNodes][1,catBFaces]
      catBFaceRightNodes=grid[BFaceNodes][2,catBFaces]
      BNodes=union(catBFaceLeftNodes,catBFaceRightNodes)
      bvals[i]=BNodes
   end
   
   return bvals
end

The user could then simply call a solution vector as indexed by the boundary face regions, e.g. if we want to extract the boundary values of a solution U associated with a system system on regions [1,3], we would simply execute

>bvals=bvalinds(system,[1,3]);
>U[1,bvals[1]] #boundary values for face region 1
>U[1,bvals[2]] #boundary values for face region 3

On the spot, I am not quite sure yet how we would extend that for arbitrary species and boundary species though and it might be worth thinking about if we should preserve some kind of ordering or orientation of the boundary nodes.

Postprocess extraction of gradients at boundaries

Hi,

I have a 2.5D (cylindrical axisymmetric) domain, over which I solve Laplace's equation for heat transfer. I have obtained the temperature gradient by doing:
∇T = nodeflux(sys, T) , which returns a matrix with the same dimensions as T.

I would like to obtain the temperature gradient on the West boundary, but am unsure what indexing yields this. I have surmised that grid[BFaceNodes], grid[BFaceEdges] or grid[BEdgeNodes] contain the indices of all boundary nodes.

Further to this, once I have obtained ∇T_west, would I be able to simply input this as a boundary condition to another PDE on the same domain? For example:
boundary_neumann!(f, p, bnode; species = 1, region = 4, value = ∇T_west)

Error with user-defined data

I wanted to pass parameters with something like

data = (D=1.0,)
phys = VoronoiFVM.Physics(flux=flux!,data=data)

but I get an error

MethodError: no method matching isdata(::NamedTuple{(:D,), Tuple{Float64}})
Closest candidates are:
  isdata(::VoronoiFVM.NoData) at /home/mjyshin/.julia/packages/VoronoiFVM/sd93r/src/vfvm_physics.jl:28
  isdata(::VoronoiFVM.AbstractData) at /home/mjyshin/.julia/packages/VoronoiFVM/sd93r/src/vfvm_physics.jl:29

Stacktrace:
 [1] VoronoiFVM.Physics(; num_species::Int64, data::NamedTuple{(:D,), Tuple{Float64}}, flux::Function, reaction::Function, storage::Function, source::Function, breaction::Function, bstorage::Function)
   @ VoronoiFVM ~/.julia/packages/VoronoiFVM/sd93r/src/vfvm_physics.jl:121
 [2] top-level scope
   @ In[77]:1
 [3] eval
   @ ./boot.jl:360 [inlined]
 [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1094

I have not gotten this error before, so I was wondering if there's something with a recent update that might have prompted it. Thank you in advance!

Could please help to solve this problem?

Could you please help to solve it in Julia with VoronoiFVM?
I can do it in ansys fluent , but I want learn Julia.

oscillating flow through pipe (2D)
it is a transient compressible, viscous 2D flow problem. also, conventional heat transfer considers from the pipe wall.
it is a convention diffusion problem. at both ends, there are porous heat exchangers.

Thank you

simplified interface

Hi,

I am really interested in using your package for solving McKean-Vlasov equations describing the distribution of mean field of particles. All examples I have in mind are like (say in 2d), x=(x1, x2)

$$\partial_tp = \partial_{x1}(F_1(x,p)p) + \partial_{x2}(F_2(x,p)p) + \Delta p$$

For example F_1(x,p) = x1^2 + sum(lambda(x), p)

Let us skip the p dependency on the vector field for now.

However, I find there is quite some boilerplate to implement in order to use your package.

1. Would it be possible to borrow the definition of the problem like in Gridap? I understand if this is aking for too much ;)

Most, if not all, your examples are nonlinear reaction diffusion equations like

du/dt = ∇(c(x,u) + d(x,u) ∇ NL(x,u)) + f(x, u)

2. Could you not define a struct asking for the functions c,d,NL,f, the discretization (like upwind or other), the boundary conditions and return a problem that implement the discretization. Even if this is a subset of what your package can handle, I'd say it would be helpful.

That you for your valuable tool

Question about evolving mesh

Hi, I am thinking of using VoronoiFVM.jl for solving a set of different PDEs governing the evolution of snow on the ground in 1D. In order to follow the ice-matrix in snow, it is usual to solve diffusion processes and water liquid percolation on a fix grid and then to update the grid (by moving nodes as snow settles) and iterate again on the new grid and so on.
I would like to know if according to you, VoronoiFVM could be applied in such a configuration.

Thank you by advance!

Replace VoronoiFVM.plot by ExtendableGrids.plot ?

VoronoiFVM.jl/examples/Example120_ThreeRegions1D.jl

Line 26 in 30f7159

VoronoiFVM.plot(Plotter,grid,p=p)

Should this be replaced with

ExtendableGrids.plot()

according to v0.8?

Unfortunately
Example120_ThreeRegions1D.main(;n=30,Plotter=Plots,plot_grid=false, verbose=false,unknown_storage=:sparse)
returns
MethodError: no method matching plot(::Module, ::ExtendableGrids.ExtendableGrid{Float64,Int32}; p=Plot{Plots.GRbackend() n=0})

for both VoronoiFVM.plot() and ExtendableGrids.plot() .

Update _initialize! to reflect new boundary condition API

As of the last API update, all boundary condition specifications can be pooled into one function and passed to the System constructor using the bcondition keyword argument.

However, the old boundary condition handling using boundary_factor and boundary_value is of course still present not only in eval_and_assemble! among

VoronoiFVM.jl/src/vfvm_solver.jl

Line 561 in 70739a1

F[ispec,K]+=boundary_factor*(U[ispec,K]-boundary_value)

but also in _initialize! e.g.

VoronoiFVM.jl/src/vfvm_system.jl

Line 931 in 70739a1

if system.boundary_factors[ispec,ibreg]≈ Dirichlet

which is invoked in every solve! call through eval_and_assemble!.

It might be sensible to commit fully to the new boundary condition handling (which pleases me very much btw, thank you for that!) because fields like boundary_factor which are assigned once within the System construcutor can lead to problems for some narrow applications where grid-dependent boundary conditions are needed.

---

Example:

I want to evaluate a convection-diffusion system on a simple simplexgrid, but set Neumann boundary conditions on one boundary face of one cell of the grid. If I want to cycle through different boundary faces in some loop, I would set up the system once where I assign a new bregion, for a 2D rectangle say ibreg=5, and specify my Neumann condition among others in some bcond function:

function bconditions(y,u,bnode)
       
        boundary_neumann!(y,u,bnode,species=1,region=5,value=1)

        # some other bconds...
        
end 

If I then make a solve! call for this system somewhere else, VoronoiFVM protests since boundary_factor will be out of bounds, having been assigned fixed when the System constructor was called.

I currently avoid this behaviour by simply setting one entry in grid[BFaceRegions] to 5 before calling the constructor and resetting it afterward.
I am sure this is not how it is intended to be used. I have not tested this further, but if the _initialize! functions were to simply invoke the breaction, being an alias for bcondition, and boundary_factor, boundary_value were made obsolete, that would at least avoid the bounds error here and make the code more consistent.

---

Apologies for the wall of text 🙃

Better documentation of test functions

See https://julialang.zulipchat.com/#narrow/stream/379007-voronoifvm.2Ejl/topic/stream.20events/near/351741229

VoronoiFVM.norm for transient solutions

Hi !

I am having a bit of trouble using the norm function for time-dependent solutions.
For example VoronoiFVM.norm(sys, u(t), 2), where u is a transient solution, gives out an error message (see below).
It might be useful for error and convergence analysis of parabolic problems.

Error Message :

MethodError: no method matching _initialize_inactive_dof!(::Vector{Float64}, ::VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}, Matrix{Float64}})

` make_all` in `docs/make.jl` ignores all input variables.

We have no control about the generation of examples and notebooks via the switches with_notebooks and with_examples at the end of the file.

couple two physics in two domains

Hi,
I am new to VoronoiFVM.jl. After one week try, I found it is really fantastic at solving diffusion problems, it is faster than similar problem I solved in c++ using FEM. So, I decide to try further.
Currently I have the time-dependent mass diffusion equation (c) in 1D domain, and heat conduction equation (T) in a 2D domain. the top boundary of the 2D domain is actually the 1D domain where mass diffusion happens. The equation are coupled (mass diffusion coefficient depends on temperature e.g., D(T), mass diffusion also generate heat, e.g., Q(c)) .

+=============mass diffusion============+
+++++++++++++++++++++++++++++++++++++++
++++++++++++++ heat conduction+++++++++++
+++++++++++++++++++++++++++++++++++++++

I first solve the 2D heat equation and get temperature T, then I solve the mass diffusion equation and get generated heat Q, Can you tell me how can I pass the T, and Q to each other in VoronoiFVM.jl?

Consolidate use of Triangle

The use of Triangle in the various Julia packages needs to be consolidated.

We have TriangleMesh by @konsim83 and Triangle by @furio . Both hide the full functionality of Triangle behind their particular mesh structure.
In this package I need yet another particular mesh structure (I need triangles, triangle attributes, segments, segment markers). I was not able to access these data with either package.
I think there should be one module e.g. TriangleRaw which gives full Julianic access to Triangle which everyone can use. Based on @konsim83 's code (which has some good architectural ideas within this respect) I made a prototype for this which sits in the subdirectory src/triangle of this package. Before making yet another Triangle package out of this I suggest to try to have some consensus.
Same IMHO would be good for TetGen and possibly other mesh generators. CC to @SimonDanisch and @ChrisRackauckas .

Jürgen (Fuhrmann)

PS: What is the appropriate place for discussing this ? It is on the intersection of geometry and PDEs...

interface to DDE problems in DifferentialEquations.jl

The interface to ODE solvers in DifferentialEquations.jl is fantastic.
What would be really great is to have an interface for DDE problems as well.
Is it feasible or too much to ask for? Thanks!

Bug in index call of quantity

VoronoiFVM.jl/src/vfvm_quantities.jl

Line 292 in bf823b3

Base.getindex(u::AbstractNodeData,q::AbstractQuantity)=@inbounds f[q[u.geom]]

It should say

u[q[u.geom]], since f is not defined.

Fix allocation in interpolation

More dots in

VoronoiFVM.jl/src/vfvm_transientsolution.jl

Line 67 in d6c8b20

retval .= a * sol[idx - 1] + b * sol[idx]

...

Density-dependent diffusion

Hello, I had asked a question a while ago about time and space dependence in the flux and reaction terms. I'm not sure if that has been resolved, but I was also wondering how one would go about modeling a flux term like j = -D(u/K)ᵐ∇u? This was my attempt (by averaging u[1,1] and u[1,2]), but I'm not sure about it:

function g!(g,u,∂ω,θ)
    ū = (u[1,1]+u[1,2])/2
    g[1] = -θ.D*(ū/θ.K)^θ.m*(u[1,2]-u[1,1])
end

Or more generally what would I do for j = -f(u)∇u? Also, since j⋅n̂ is approximated by g(uₖ,uₗ)/|xₖ-xₗ|, wouldn't we need to multiply the u by something like (abs ∘ -)(∂ω.coord[1:2]...)? Thank you in advance!