The source in Poisson's equation determines the curvature of the solution. Neumann boundary conditions determine the slope at the boundaries. There are still infinitely many solutions that satisfy those conditions.
> On Jul 23, 2016, at 7:47 AM, Gopalakrishnan, Krishnakumar > <[email protected]> wrote: > > 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: [email protected] [mailto:[email protected]] On Behalf Of > Guyer, Jonathan E. Dr. (Fed) > Sent: Thursday, July 14, 2016 5:26 PM > To: FIPY <[email protected]> > 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) >> <[email protected]> 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 >>> <[email protected]> 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 >>> [email protected] >>> http://www.ctcms.nist.gov/fipy >>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >> >> >> _______________________________________________ >> fipy mailing list >> [email protected] >> http://www.ctcms.nist.gov/fipy >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
