usnistgov / fipy

FiPy is a Finite Volume PDE solver written in Python
http://pages.nist.gov/fipy/en/latest
Other
489 stars 148 forks source link

add a baffle or stone to the simulation #978

Closed QidiMa closed 7 months ago

QidiMa commented 7 months ago

Hi all,

It should be quite obvious but I'm a beginner with Fipy, apologies! Here is my simulation code, and I want to add a baffle or stone to the simulation .What should I do.

Thanks a lot! Best regards,

###
from fipy import CellVariable, Grid2D, TransientTerm, DiffusionTerm, ExponentialConvectionTerm, Viewer
import fipy as fp
import numpy as np

# create mesh
nx, ny = 100, 100
dx, dy = 1., 1.
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

# velocity
u = np.ones((100, 100))*10
v = np.ones((100, 100))
u_var = CellVariable(mesh=mesh, value=u.flatten())  # 10000
v_var = CellVariable(mesh=mesh, value=v.flatten())  # 10000
u_face_var = u_var.arithmeticFaceValue  # 20200
v_face_var = v_var.arithmeticFaceValue  # 20200
velocity = fp.FaceVariable(mesh=mesh, value=[u_face_var, v_face_var], rank=1)

# define concentration
c = CellVariable(name="Concentration", mesh=mesh, value=0.)

# boundary condition
c.constrain(0., mesh.exteriorFaces)

# parameters
# diffusion coefficient
D = 1

# time setting
dt = 0.1
steps = 50

# source setting
source_duration = 10
source_value = 10.0
r = 4.0
center_x, center_y = nx / 2, ny / 2

# viewer setting
viewer = Viewer(vars=c)

# time iteration
for step in range(steps):
    if step < source_duration:
        source_term = source_value * ((mesh.x - center_x) ** 2 + (mesh.y - center_y) ** 2 <= r ** 2)
    else:
        source_term = 0.

    # define equation
    eq = (TransientTerm(var=c) == DiffusionTerm(coeff=D, var=c) - ExponentialConvectionTerm(coeff=velocity, var=c) + source_term)
    eq.solve(var=c, dt=dt)
    # save picture
    if step % 3 == 0:
        viewer.plot()
    # viewer.plot(filename=f"{step}.png")