In your initial condition, w is zero everywhere, correct? In that case, your PDE is indistinguishable from Laplace's equation and so has infinitely many solutions.
You might try setting your initial condition of w to something nonzero to see if that helps initial convergence. > On Jul 26, 2016, at 10:10 AM, Gopalakrishnan, Krishnakumar > <krishnaku...@imperial.ac.uk> wrote: > > Hi Jon, > > Thanks a lot for replying. > > Based a reply to one of our earlier posts on this mailing list on 26 May 2016 > *, we understand that: > > If the right-hand-side (i.e. source-term) of the elliptic PDE is a function > of the field-variable, then this uniquely pins down the solution vector > (even though we have only a Neumann BC). I am quoting the relevant sentence > from that post below. > > "However, the form of the source is important here because if the integral of > the source is non-zero, then there is no steady solution. If the integral > _is_ zero and the source is not a function of \phi, then the system admits an > infinite family of solutions, all shifted by a constant." - Raymond Smith > > * http://thread.gmane.org/gmane.comp.python.fipy/4110/focus=4111 > > > > We tried this approach on a standalone elliptic PDE with implicit source term > (albeit without the coupling into the 2D domain with anisotropic diffusion), > and it was successful. However, when we modify your gist code, it throws up > errors with Zero-flux Neumann BC on the left boundary (even though the > f(w,a) is now replaced with ImplicitSourceTerm (coeff=w**3.0). I assume that > this should have uniquely constrained the solution, as it did in our > standalone implementation. Kindly correct us if our understanding is lacking > in this context. > > We definitely need to implement the zero-flux Neumann BC anyway, with the > crazy coupling described in the original post. Can you recommend a solution > for implementing this ? > > > Best Regards, > > Ian & Krishna > > > -----Original Message----- > From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of > Guyer, Jonathan E. Dr. (Fed) > Sent: Tuesday, July 26, 2016 1:16 PM > To: FIPY <FIPY@nist.gov> > Subject: Re: Coupling Boundary Condition of one PDE with source term of > another PDE > > 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 >> <krishnaku...@imperial.ac.uk> 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: 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 ] > > > _______________________________________________ > 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 ]
