[note: I sent this to Mohammad more than a week ago and only learned now that
it didn't go to the list]
Mohammad -
A couple of things:
- dc/dn = 0 is *not* a no flux condition. This is a really common assumption,
and it's really wrong when there's strong convection. Your flux is $U_f \phi -
D \nabla \phi$.
- FiPy is a cell-centered finite volume solver. It's solving the divergence
theorem on every cell. As a result, no flux boundary conditions are the
default. Not only do you not need to set a constraint, any constraint you set
will probably cause the boundaries to *not* be no flux.
Resolving these issues don't fix your problem, though. Ultimately, I think
you're seeing phi exceed 200 because Uf is not divergence-free. You can see
this by declaring a viewer
>>> viewer = fp.Viewer(vars=Uf.divergence, datamin=-0.0001, datamax=0.0001)
The divergence at the outlet river is not a problem, but everywhere else, it's
saying that your lake water is compressible.
As far as I understand your script and data file, this is happening because
FakeLake_SED_WU_Test.nc contains a 3D flow field and you're approximating it as
a 2D field by averaging over the six depth layers. Even so averaged, there are
still vertical flows, which amount to sources and sinks when you project it in
2D. These sources and sinks seem particularly strong around the neck of "top
river", which is where phi is going high.
Frankly, I would just model the 3D flow; it's only 6x the cells you're modeling
right now. Hopefully the full 3D flow data set is divergence free.
Other than that:
- I would only constrain phi.faceGrad to zero at the outlet river
- Ideally, your inlet and outlet masks should be defined parametrically or,
even better, be explicitly identified by your mesh generator. "mask2 = (YY <
4692200)" grabs a little bit of the bottom left river as well as some of the
"bank" of the bottom right river. I'm not positive, but I think it doesn't
completely capture the inlet of the bottom right river. I don't know what SMS
provides, but you should try to make it tell you where the boundaries are,
rather than hard wiring things like "- 50" and "< 4692200".
- Don't redefine the constraints and equation at every time step. It definitely
slows things down, and probably leaks memory, too.
- I'm not sure if it matters, but the normal component of Uf should be zero on
the exterior boundaries, but not at the river inlet/outlets.
- Jon
> On Nov 24, 2019, at 3:59 PM, Mohammad Madani <mada...@uwindsor.ca> wrote:
>
> Thank you for helping me with the mesh.
> Please find attached files. I created the mesh and set up the model.
> Regarding the boundary condition I have three conditions.
>
> 1 - no flux in all exterior faces (dc/dn = 0) except in input rivers.
> 2 - constant concentration of 200 in top river
> 3- constant concentration of 100 in small river in bottom right.
>
> I created those boundary conditions like this:
> XX, YY = mesh.faceCenters
> mask = YY > np.max(YY) - 50 # top river
> mask2 = YY < 4692200 # bottom right river
> phi.faceGrad.constrain(0., where=mesh.exteriorFaces) # no flux condition
> phi.constrain(200., where=mesh.exteriorFaces & mask)
> phi.constrain(100., where=mesh.exteriorFaces & mask2)
>
> Are they correct?
> I notice there is a problem. With this inputs I should not have concentration
> above 200 in my domain however I got values around 250-280 ?
>
> I somehow found out what the problem maybe is but I cannot fix it. It is
> related to the convection coefficient or field velocities. if I use constant
> (-1 , -1 ) in line 156 !
>
> U_t = -1* np.ones(shape=(2, mesh.numberOfFaces))
> Uf.setValue(U_t)
> the model will run without getting any value higher than 200 but when I use
> the actual Ux and Uy component of velocities, it gives me unrealistic values.
> Could you please help me to fix this.
>
> Attached are my code and one nc file that contains mesh node and element
> information, temperature, and velocity components.
>
>
> Thanks
> Mohammad
>
> ------ Original Message ------
> From: "Guyer, Jonathan E. Dr. (Fed) via fipy" <fipy@nist.gov>
> To: "FIPY" <FIPY@nist.gov>
> Sent: 2019-11-22 11:30:19 AM
> Subject: Re: Need help to set up the model
>
>> Mohammad -
>>
>> Welcome to FiPy.
>>
>> It's not a problem to mix triangles and quadrilaterals. FiPy is designed for
>> this. The complication is that FiPy constructs cells from faces and faces
>> are defined by vertices (nodes), so we need to reconstruct the face
>> definitions from your cell_node_id.
>>
>> Cribbing from gmshMesh.py:
>> ============
>> import fipy as fp
>> from fipy.meshes.mesh2D import Mesh2D
>> from fipy.tools import numerix as nx
>>
>> node_coordinates = fp.numerix.array(
>> ((0., 0.),
>> (1., 0.),
>> (0., 1.),
>> (1., 1.),
>> (2., 0.),
>> (3., 0.),
>> (2., 1.),
>> (3., 1.),
>> (2., 2.),
>> (3., 2.)
>> )
>> )
>>
>> cell_node_id = fp.numerix.MA.masked_less((
>> (0, 1, 2, -1),
>> (1, 3, 2, -1),
>> (1, 4, 6, 3),
>> (4, 5, 7, -1),
>> (4, 7, 6, -1),
>> (6, 7, 9, 8)
>> ), value=0.)
>>
>> currNumFaces = 0
>> cellsToFaces = nx.ones(cell_node_id.shape, long) * -1
>> facesDict = {}
>> uniqueFaces = []
>>
>> for cellIdx, cell in enumerate(cell_node_id):
>> cell = cell.compressed()
>> faces = fp.numerix.concatenate((cell[..., nx.newaxis],
>> fp.numerix.roll(cell, -1)[..., nx.newaxis]),
>> axis=1)
>> for faceIdx, face in enumerate(faces):
>> keyStr = ' '.join([str(x) for x in sorted(face)])
>>
>> if keyStr in facesDict:
>> cellsToFaces[cellIdx][faceIdx] = facesDict[keyStr]
>> else: # new face
>> facesDict[keyStr] = currNumFaces
>> cellsToFaces[cellIdx][faceIdx] = currNumFaces
>> uniqueFaces.append(list(face))
>> currNumFaces += 1
>>
>> # pad short faces with -1
>> maxFaceLen = max([len(f) for f in uniqueFaces])
>> uniqueFaces = [[-1] * (maxFaceLen - len(f)) + f for f in uniqueFaces]
>>
>> facesToVertices = nx.array(uniqueFaces, dtype=int)
>>
>> mesh = Mesh2D(vertexCoords=node_coordinates.swapaxes(0, 1).copy('C'),
>> faceVertexIDs=facesToVertices.swapaxes(0, 1)[::-1],
>> cellFaceIDs=cellsToFaces.swapaxes(0, 1).copy('C'))
>> ============
>>
>>
>> In order to help with the boundary conditions, I need to know more about how
>> define them, both their location and their values.
>>
>> - Jon
>>
>>> On Nov 21, 2019, at 10:30 PM, Mohammad Madani <mada...@uwindsor.ca> wrote:
>>>
>>> Hello Everyone,
>>> I sent below email to Dr. Guyer and he advised me to join the mailing list
>>> and bring my questions here. Thank you!
>>>
>>> Attached is a python file that I tried to modify one of the examples based
>>> on my need.
>>> First of all, I have a big problem in creating my mesh. I could not find
>>> how to import my mesh to Fipy from the manual. Would you please help me
>>> with that. The one showed in the attached file is just a template and it is
>>> not my mesh. For my mesh, I have node_coordinates and cell_node_id: The
>>> mesh was created using SMS software.
>>> node_coordinates:
>>> array([[ 360980. , 4709150. ],
>>> [ 361720. , 4709170. ],
>>> [ 360610. , 4709140. ],
>>> ...,
>>> [ 354182. , 4692112. ],
>>> [ 353605. , 4692098.2],
>>> [ 353914. , 4691819. ]])
>>>
>>> cell_node_id:
>>> array([[ 0, 3, 4, --],
>>> [ 2, 7, 3, 0],
>>> [ 5, 1, 4, 9],
>>> ...,
>>> [1005, 1009, 1011, 1008],
>>> [1006, 1010, 1012, 1009],
>>> [1013, 1011, 1009, 1012]])
>>>
>>> my mesh has both triangles and quadrilateral cells. if that makes problem I
>>> can create my mesh just with triangles
>>>
>>> After creating my mesh, I know how to set up the model as I did in the
>>> attached file, but still I don't know how to set up Boundary Condition.
>>> from node_string 1 and 2 (below figure), variable concentration over time
>>> introduced to the lake.
>>>
>>> Any comments and help would be appreciated.
>>>
>>> Thanks
>>> Mohammad
>>>
>>> ------ Forwarded Message ------
>>> From: "Mohammad Madani" <mada...@uwindsor.ca>
>>> To: jonathan.gu...@nist.gov
>>> Sent: 2019-11-14 11:50:31 AM
>>> Subject: Need help to set up the model
>>>
>>> Hello Dr. Guyer
>>> I am PhD student from University of Windsor, Canada. I am working on
>>> microbial modelling in the lake and would like to use Fipy to solve the E.
>>> coli transport equation. the governing equation is
>>>
>>> <wnxflldd.png>
>>>
>>> where C is the E. coli concentration and k is the decay rate. k is function
>>> of temperature and solar radiation which I have the time series of solar
>>> radiation. I also have water temperature for each cell from the
>>> hydrodynamic model.
>>> lets say the k value can be calculate as:
>>>
>>> k = k0 * 10 ^( 20 - T) * Solar_Radiation
>>>
>>> where k0 is constant and T is temperature.
>>>
>>> I also have time series of velocities (u, v and w) from hydrodynamic model
>>> for each cell. Here is my mesh:
>>> <rjrrcjsy.png>
>>> which node-string 1 and 2 are input and 3 is the output of the lake. E coli
>>> concentration from the inputs can be provided as time series too.
>>> Can I solve this problem with Fipy? would you please help me to set up the
>>> python code to do it.
>>>
>>> looking forward to hear from you soon.
>>> Thank you,
>>>
>>> <image.jpg>
>>> Mohammad Madani
>>> PhD Candidate
>>> Department of Civil & Environmental Eng.
>>> CEI, Room 3084, Desk-I5
>>> 401 Sunset Ave. Windsor ON N9B 3P4
>>> Mobile: 519-992-8267
>>> Email: mada...@uwindsor.ca
>>> uwindsor.ca/engineering
>>> twitter.com/UWindsorENG
>>> facebook.com/UWindsorEngineering
>>>
>>> <Fipy_Lake_Model.py>_______________________________________________
>>> fipy mailing list
>>> fipy@nist.gov
>>> http://www.ctcms.nist.gov/fipy
>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>
>>
>> _______________________________________________
>> fipy mailing list
>> fipy@nist.gov
>> http://www.ctcms.nist.gov/fipy
>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> <Fipy_Lake_Model_V1.py><FakeLake_SED_WU_Test.nc>
from fipy import FaceVariable, CellVariable, TransientTerm, DiffusionTerm, PowerLawConvectionTerm
from fipy import ImplicitSourceTerm, MatplotlibViewer, DefaultAsymmetricSolver, Variable, ExponentialConvectionTerm
from fipy.solvers import LinearPCGSolver, LinearLUSolver, GeneralSolver
import fipy as fp
from fipy.meshes.mesh2D import Mesh2D
from fipy.tools import numerix as nx
from netCDF4 import Dataset
import numpy as np
from matplotlib import pylab
nc_file = 'FakeLake_SED_WU_Test.nc'
f = Dataset(nc_file)
cell_node = f.variables['cell_node']
cell_X = f.variables['cell_X']
cell_Y = f.variables['cell_Y']
cell_Z = f.variables['cell_Zb']
cell_Nvert = f.variables['cell_Nvert']
node_X = f.variables['node_X']
node_Y = f.variables['node_Y']
node_Z = f.variables['node_Zb']
num_node = node_X[:].shape[0]
layerface_Z = f.variables['layerface_Z']
Ux = f.variables['V_x']
Uy = f.variables['V_y']
Uz = f.variables['W']
Time = f.variables['ResTime']
# convert simulation time to second
if Time.units == 'hours':
Time = Time[:] * 3600
def interpolate_data(data, t, Time):
# the first index of data should be time !
dt = Time - t
indx = np.where(dt == 0)[0]
if not indx.size == 0:
data_interpolated = data[indx, ...]
else:
first = np.where(dt < 0)[0][-1]
last = np.where(dt > 0)[0][0]
ratio = (t - Time[first]) / (Time[last] - Time[first])
data_interpolated = data[last, ...] * ratio + data[first, ...] * (1 - ratio)
return data_interpolated
U = np.zeros((13, 1303, 6, 2))
U[:, :, :, 0] = Ux[:].reshape(13, 1303, 6)
U[:, :, :, 1] = Uy[:].reshape(13, 1303, 6)
# depth average values
U = np.nanmean(U, axis=2)
TEMP = np.nanmean(f.variables['TEMP'][:].reshape(13, 1303, 6), axis=2)
node_coordinates = np.vstack((node_X[:], node_Y[:])).T
cell_node_id = fp.numerix.MA.masked_less(cell_node[:] - 1, value=0.)
currNumFaces = 0
# long = 5
cellsToFaces = nx.ones(cell_node_id.shape, long) * -1
facesDict = {}
uniqueFaces = []
for cellIdx, cell in enumerate(cell_node_id):
cell = cell.compressed()
faces = fp.numerix.concatenate((cell[..., nx.newaxis],
fp.numerix.roll(cell, -1)[..., nx.newaxis]),
axis=1)
for faceIdx, face in enumerate(faces):
keyStr = ' '.join([str(x) for x in sorted(face)])
if keyStr in facesDict:
cellsToFaces[cellIdx][faceIdx] = facesDict[keyStr]
else: # new face
facesDict[keyStr] = currNumFaces
cellsToFaces[cellIdx][faceIdx] = currNumFaces
uniqueFaces.append(list(face))
currNumFaces += 1
# pad short faces with -1
maxFaceLen = max([len(f) for f in uniqueFaces])
uniqueFaces = [[-1] * (maxFaceLen - len(f)) + f for f in uniqueFaces]
facesToVertices = nx.array(uniqueFaces, dtype=int)
# Define mesh
mesh = Mesh2D(vertexCoords=node_coordinates.swapaxes(0, 1).copy('C'),
faceVertexIDs=facesToVertices.swapaxes(0, 1)[::-1],
cellFaceIDs=cellsToFaces.swapaxes(0, 1).copy('C'))
# Define variables
phi = CellVariable(name=r"$E.coli$ Concentration",
mesh=mesh,
value=0.,
hasOld=True)
K = CellVariable(name=r'Decay rate ($d^{-1}$)', mesh=mesh, value=0.)
Uf = FaceVariable(mesh=mesh, rank=1)
Uc = CellVariable(mesh=mesh, rank=1)
Ux = CellVariable(mesh=mesh, value=0.)
Uy = CellVariable(mesh=mesh, value=0.)
Tmp = CellVariable(mesh=mesh, name='Temperature', value=0.)
time = Variable()
D = 0.001
viewer = None
if __name__ == '__main__':
try:
# viewer = fp.MatplotlibStreamViewer(vars=Uf)
viewer = fp.Viewer(vars=Uf.divergence, title=r"$\nabla\cdot\vec{U}$", datamin=-0.0001, datamax=0.0001)
# viewer = MatplotlibViewer(vars=phi,
# title=r'$E. coli$ Concentration (CFU/100 mL)',
# datamin=0.,
# datamax=300.,
# colorbar='horizontal',
# figaspect=0.7,
# cmap=pylab.cm.jet) # pylab.cm.OrRd)
# viewer2 = MatplotlibViewer(vars=Ux,
# title=r'Decay rate ($d^-1$)',
# datamin=-1.,
# datamax=1.,
# colorbar='horizontal',
# figaspect=0.7,
# cmap=pylab.cm.jet) # pylab.cm.OrRd)
# viewer.plotMesh(mesh)
except:
print("Unable to create a viewer for an irregular mesh (try Matplotlib2DViewer or MayaviViewer)")
phi.setValue(100.)
timeStepDuration = 360. # second
XX, YY = mesh.faceCenters
mask = YY > np.max(YY) - 50 # top river
mask2 = (YY < 4692200) & (XX > 365000) # bottom right river
mask3 = (YY < 4693500) & (XX < 354000) # bottom left river
elapsedTime = 0
# Uf.constrain(0., where=mesh.exteriorFaces & ~(mask | mask2 | mask3))
faces = fp.FaceVariable(mesh=mesh, value=1. * (mesh.exteriorFaces & ~(mask | mask2 | mask3)))
fp.Viewer(vars=(faces.divergence != 0.) * 1., title="Boundary Cells").plot()
# phi.faceGrad.constrain(0., where=mesh.exteriorFaces & mask3) # outflow condition
phi.faceGrad.constrain(0., where=mesh.exteriorFaces) # outflow condition
phi.constrain(200., where=mesh.exteriorFaces & mask)
phi.constrain(100., where=mesh.exteriorFaces & mask2)
# phi.faceGrad.constrain(0., mesh.exteriorFaces & ~mask)
eq = TransientTerm(coeff=1.) + PowerLawConvectionTerm(coeff=Uf) == \
DiffusionTerm(coeff=D) + ImplicitSourceTerm(coeff=K)
while time.value < Time[-1]:
print("Time", '*' * (60 - len(str(time))), time, "(s)")
# if hasattr(phi, 'faceConstraints'):
# del phi.faceConstraints
Tmp.setValue(interpolate_data(TEMP, time.value, Time))
K.setValue(-0.001 * 10 ** (20 - Tmp))
Ux.setValue(interpolate_data(U[:, :, 0], time.value, Time))
Uy.setValue(interpolate_data(U[:, :, 1], time.value, Time))
# UU = CellVariable(mesh=mesh, value=U[1,:,1,:].T)
# UU = FaceVariable(mesh=mesh, value=np.mean(U[step, mesh.faceCellIDs, 1, :], axis=0).T)
# calculate velocity in the cell faces
U_t = np.vstack((Ux.arithmeticFaceValue, Uy.arithmeticFaceValue))
# U_t = -1* np.ones(shape=(2, mesh.numberOfFaces))
Uf.setValue(U_t)
Uc[0] = Ux
Uc[1] = Uy
residual = 1.0
phi.updateOld()
sweep_iteration = 0
while residual > 1e-5:
sweep_iteration += 1
residual = eq.sweep(var=phi,
dt=timeStepDuration,
solver=DefaultAsymmetricSolver(tolerance=1.e-15, iterations=10000))
if sweep_iteration > 100:
break
time.setValue(time() + timeStepDuration)
print('Maximum Concentration = {0:.2f} CFU/100mL'.format(phi.value.max()))
if viewer is not None:
viewer.plot()
# viewer2.plot()
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
