Hi Jon, Thanks a lot for your valuable code in getting us onto the right track.
However, there is an immediate problem when we start to modify the code slightly. If the BC1 of PDE1 is changed from Dirichelet to (zero-flux) Neumann, the solver stagnates. In addition, we get the error 'ScalarQuantityOutOfRangeWarning: A scalar quantity became too small or too large to continue computing. Iterations: 15. Relative error: 1.63296e+16' when the interpreter gets to the line, res1 = eq1.sweep(a). The results have a lot of NaNs and a few other error and warning messages follow. This is the case in both the crazy_coupling1 and crazy_coupling2 codes. Zero-flux Neumann is a valid BC for PDE1 since the PDE is implicit, i.e. the source term on the right is defined by f(w, a), where a is the CellVariable. This should pin down the value at the left boundary and hence make the physics valid. How can we massage your gist code to be amenable to this? Best Regards, Krishna -----Original Message----- From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Guyer, Jonathan E. Dr. (Fed) Sent: Thursday, July 14, 2016 5:26 PM To: FIPY <FIPY@nist.gov> Subject: Re: Coupling Boundary Condition of one PDE with source term of another PDE My gut reaction is that the second implementation is probably better, as it should be possible to make things more implicit and coupled, even though you're "wasting" calculation of PDE1 over most of the 2D domain where you don't care about it. > On Jul 14, 2016, at 11:18 AM, Guyer, Jonathan E. Dr. (Fed) > <jonathan.gu...@nist.gov> wrote: > >> What’s the best way to implement this problem in FiPy? > > Don't! > > > Assuming you won't take that advise, I've posted a couple of attempts at this > problem at: > > https://gist.github.com/guyer/bb199559c00f6047d466daa18554d83d > > > >> On Jul 9, 2016, at 1:37 PM, Gopalakrishnan, Krishnakumar >> <k.gopalakrishna...@imperial.ac.uk> wrote: >> >> Hi, >> >> **We’re sorry, the previous posting omitted the complete definition of the >> boundary condition for the 2nd PDE, and also a couple of typos are fixed >> below. >> >> We have a system of equations, wherein the BC of one PDE couples with the >> source term of another PDE. >> >> We have a regular 2D unit grid in x and y. >> >> There are two PDEs to be solved: >> a) The first PDE (elliptic diffusion problem) is defined only at y = 1, >> acting along the x-axis (i.e. it acts in the x-direction and only along the >> top of the cartesian grid). This x-axis is discretized with a fixed >> grid-spacing, generating a finite number of nodes. Let this set of >> nodes/co-ordinates be represented by ‘X’. >> b) The second PDE is a time-varying diffusion problem. This is defined >> only along the y-axis, but for all x-nodes (i.e. for all ‘X’), where the 1st >> PDE is being solved. >> >> Mathematically (plain-text version) : >> >> PDE1: divergence(S * grad(a)) = f(w,a) >> BC1: a = 0 , at x = 0 (dirichelet) >> BC2: d_a/dx = 1 at x = 1 (Neumann) >> >> D = faceVariable(rank=2, value = 1.0) # Rank 2 tensor for anisotropic >> diffusion (i.e. allowed only along the y-axis) >> >> PDE2: partial_B/partial_t = divergence( [[[0, 0], [0, D]]] * grad(B) ) >> BC1: B = 0, at y = 0 at all ‘X’ (dirichelet), i.e. along the bottom face >> BC2: d_B/dy = g(w,a) (Neumann) at y =1, i.e. along the top face of the grid >> >> f and g are linear functions in ‘w(x,y,t)’ and ‘a(x,y,t)’. B is defined in >> the 2D grid as ‘B(x,y,t)’. >> >> Importantly, w = B at y = 1, for all ‘X’; i.e., the BC2 of the 2nd PDE >> couples with the Implicit Source Term of the 1st PDE along the top face of >> the Cartesian mesh. >> >> What’s the best way to implement this problem in FiPy? >> >> Best Regards, >> >> Krishna >> >> >> PS: Since the problem is mathematically harder to express in plain-text, >> here is a HTML formatted version of the PDEs and BCs. >> PDE1: <image002.png> >> BC1: <image004.png> (Dirichelet) >> BC2: <image006.png> (Neumann) >> >> D = faceVariable(mesh=pos_p2d_mesh, rank=2, value = 1.0) # Rank 2 >> tensor for anisotropic diffusion (i.e. allowed only along the y-axis) >> >> PDE2: <image008.png> = <image010.png> >> BC1 is: <image012.png> i.e. along the bottom face of the grid >> BC2 is: <image014.png> (Neumann) at y =1, i.e. along the top face of the >> grid >> >> f and g are linear functions in <image015.png> and<image016.png>. >> <image018.png> is defined in the 2D grid as <image035.png> >> >> Importantly, <image036.png> i.e., the BC2 of the 2nd PDE couples with the >> Implicit Source Term of the 1st PDE along the top face of the Cartesian mesh. >> >> <image017.png><image024.png><image026.png><image023.png><image025.png><image022.png><image011.png><image019.png><image020.png>_______________________________________________ >> 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 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 ]
