Commit 658a4be1 authored by Andreas Dedner's avatar Andreas Dedner

added a simpler fv example

parent 63b0d7ef
import numpy as np
from mpi4py import MPI
from dune.common import FieldVector
import dune.grid as grid
from dune.grid import yaspGrid, gridFunction
# return a vector with the initial concentration by evaluating
# the grid function c0 at the barycenter of each element using
# global coordinates
def initialize(gridView, c0):
# mapper with one value on each element
mapper = gridView.mapper(lambda gt: gt.dim == gridView.dimension)
# numpy vector with data on each element
c = np.zeros(len(mapper))
for e in gridView.allPartition.elements:
center = e.geometry.center # global coordinate of cell center
c[mapper.index(e)] = c0(center) # initialize cell concentration
return c, mapper
# main evolution method - the concentration is update in place
# using 'b' at the boundary and 'u' as transport direction
# c_T^{n+1} -= dt sum_j |e^T_j| g^T_j(c_T,c_{N_j})
# where g^T_j is an upwind flux over the edges (e^T_j)_j of the element T
# using the neighboring concentration value c_{N_j}. The numerical flux is
# consistent with the linear analytical flux f(c) = uc with a given
# constant transport direction u. The time step size dt is computed from
# the maximum wave speed.
def evolve(gridView, mapper, concentration, t, b, u):
# update vector and time step
update = np.zeros(len(concentration))
dt = 1e100
# iterate over the elements
for entity in gridView.allPartition.elements:
volume = entity.geometry.volume
indexi = mapper.index(entity)
sumfactor = 0.0
# iterate over the intersections
for intersection in gridView.intersections(entity):
igeo = intersection.geometry
faceglobal = igeo.center
normalvelocity = (u(faceglobal, t) * intersection.centerUnitOuterNormal) * igeo.volume
factor = normalvelocity / volume
if factor > 0:
sumfactor += factor
outside = intersection.outside
if outside is not None: # inner intersection so use data on neighboring element
indexj = mapper.index(outside)
if indexi < indexj:
nbfactor = normalvelocity / outside.geometry.volume
# upwind concentration
c = concentration[indexj if factor < 0 else indexi]
update[indexi] -= c * factor
update[indexj] += c * nbfactor
elif intersection.boundary: # boundary intersection
c = b(faceglobal, t) if factor < 0 else concentration[indexi]
update[indexi] -= c * factor
dt = min(dt, 1.0 / sumfactor)
# reduce time step a bit
dt = 0.9 * dt
concentration += update*dt # concentration vector updated by side effect
return dt
if __name__ == "__main__":
domain = grid.cartesianDomain([0, 0], [1, 1], [40, 40])
gridView = grid.yaspGrid(domain)
@gridFunction(gridView)
def c0(x):
if x.two_norm > 0.125 and x.two_norm < 0.5:
return 1.0
else:
return 0.0
def u(x, t):
return FieldVector([1.0, 1.0])
def b(x, t):
return c0(x-t*u(x,t))
gridView.hierarchicalGrid.globalRefine(1)
concentration, mapper = initialize(gridView, c0)
vtkwriter = gridView.vtkWriter()
# make a grid function from the piecewise constant data for output
@gridFunction(gridView)
def solution(element,point):
return [concentration[mapper.index(element)]]
solution.addToVTKWriter("concentration", vtkwriter, grid.DataType.CellData)
# solution.plot() # uncomment this for matplot output
counter = 0
t, dt = 0.0, 0.0
k = 0
saveInterval = 0.1
saveStep = 0.1
# output the initial data to a vtk file
vtkwriter.write("concentration", counter, grid.OutputType.appendedraw)
counter += 1
# time loop
while t < 1.:
k += 1
dt = evolve(gridView, mapper, concentration, t, b, u)
print("k=" + str(k) + " ; dt=" + str(dt) + " ; time=" + str(t))
t += dt
if t >= saveStep:
vtkwriter.write("concentration", counter, grid.OutputType.appendedraw)
saveStep += saveInterval
counter += 1
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment