### Re: Need help to set up the model

[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 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" > To: "FIPY" > 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 impor

### 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 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" > 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 > > > > 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: > > 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, > > > Mohammad Madani > PhD Candidate > Department of Civil & Environmental Eng. > CEI, Room 3084,