Last active
June 28, 2023 08:54
-
-
Save wd15/28f1b13f9f570fdcc21b9cb22c48bf12 to your computer and use it in GitHub Desktop.
Answer to stackoverflow question, https://stackoverflow.com/questions/58695069/fipy-convection-with-a-given-velocity-field
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from fipy import Grid2D, CellVariable, FaceVariable | |
from fipy import ConvectionTerm, DiffusionTerm, TransientTerm | |
import numpy as np | |
from scipy import interpolate | |
from fipy import LinearGMRESSolver | |
dx = 0.5 | |
nx = 7 | |
dy = 0.5 | |
ny = 7 | |
xy = np.zeros((nx, ny, 2)) | |
xy[:, :, 0] = np.mgrid[0:nx, 0:ny][0] * dx | |
xy[:, :, 1] = np.mgrid[0:nx, 0:ny][1] * dy | |
u = xy[..., 0] + xy[..., 1] | |
v = xy[..., 0] * xy[..., 1] | |
print(xy.shape) | |
print(u.shape) | |
print(v.shape) | |
m = Grid2D(nx=3, ny=3, dx=1.0, dy=1.0) | |
xy_interp = np.array(m.faceCenters).swapaxes(0, 1) | |
u_interp = interpolate.griddata(xy.reshape(-1, 2), u.flatten(), xy_interp, method='cubic') | |
v_interp = interpolate.griddata(xy.reshape(-1, 2), v.flatten(), xy_interp, method='cubic') | |
var = CellVariable(mesh=m) | |
diffusion = 1.0 | |
velocity = FaceVariable(mesh=m, rank=1) | |
print(u_interp) | |
print(v_interp) | |
velocity[0, :] = u_interp | |
velocity[1, :] = v_interp | |
eqn = TransientTerm() + ConvectionTerm(velocity) == DiffusionTerm(diffusion) | |
eqn.solve(var, dt=1.0) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Thanks for the example, but the solution I get is just zeros..