Interpolated on what ? After the loop on the steps, the variable phi is still accessible, you should be able to manipulate it as a numpy array and store it or do whatever you like with it.
Thibault Bridel-Bertomeu 2017-05-18 14:37 GMT+02:00 Sergio Manzetti <sergio.manze...@fjordforsk.no>: > Thibault, the solution plotted by the Viewer looks like a decaying > exponential function. Is there any way to get out from this plot an > interpolated function, which is the result to the solved PDE? > > Thanks! > > > Sergio Manzetti > > <http://www.fjordforsk.no/logo_hr2.jpg> > > Fjordforsk AS <http://www.fjordforsk.no> > <http://www.fjordforsk.no> > Midtun > 6894 Vangsnes > Norge > Org.nr. 911 659 654 > Tlf: +47 57695621 <+47%2057%2069%2056%2021> > Økolab <http://www.oekolab.com> | Nanofactory <http://www.nanofact.no> > | AQ-Lab <http://www.aq-lab.no> | FAP <http://www.phap.no> > > > ------------------------------ > *From: *"Thibault Bridel-Bertomeu" <thibault.bridellel...@gmail.com> > *To: *"fipy" <fipy@nist.gov> > *Sent: *Thursday, May 18, 2017 2:29:23 PM > > *Subject: *Re: Problems running a simple PDE > > I invite you to take a look at : http://www.ctcms.nist.gov/ > fipy/examples/README.html for an overview of all the examples written for > fipy > and http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html > for a reminder on how the finite volume method works. > > As for your equation, if you are in 1D, you may consider it as follow : > > du/dx : is a convection term with a unit scalar coefficient, i.e. > <SpecificConvectionTerm>(coeff=(1.,), var=u) - you will find a list of > all available convection schemes here : http://www.ctcms.nist.gov/ > fipy/documentation/numerical/discret.html#convection-term > > du/dt : is a transient term that you can treat as before, i.e. > TransientTerm(var=u) > > and u^2, as a nonlinear term, one possibility is to update it during the > iterations. > > All in all, you would write for instance : > > #!/usr/bin/env python > > from fipy import * > from fipy import numerix > > nx = 50 > dx = 1. / float(nx) > > mesh = Grid1D(nx=nx,dx=dx) > X = mesh.cellCenters[0] > > phi = CellVariable(mesh=mesh, name="solution") > phi.setValue(0.5-0.5*numerix.tanh(10*(X-0.5))) > > vi = Viewer(vars=phi,datamin=0.0, datamax=1.0) > vi.plot() > > raw_input("Initialization ...") > > phi.constrain(1., mesh.facesLeft) > phi.constrain(0., mesh.facesRight) > > phi_sq = CellVariable(mesh=mesh) > phi_sq.setValue( phi*phi ) > > eq = TransientTerm(coeff=1., var=phi) + ExponentialConvectionTerm(coeff=(1.,), > var=phi) + phi_sq == 0.0 > > dt = 0.01 > steps = 100 > for step in range(steps): > eq.sweep(dt=dt) > # > phi_sq.setValue( phi * phi ) > # > vi.plot() > > raw_input("Press <return> ...") > > The sweep method allows to advance the equation of one timestep only. Then > I can update phi_sq which the next sweep will use to solve the equation. > And so on .. > > Hope this helps to understand. I can however only advise you to browse > through all the examples, you will definitely find something similar to > your case you can start from ! > > Best > > T. Bridel-Bertomeu > > > > 2017-05-18 14:10 GMT+02:00 Sergio Manzetti <sergio.manze...@fjordforsk.no> > : > >> OK, this really helped. If I wanted to change this into: >> >> >> du/dx + i*du/dt + u^2 = 0 >> >> in the following format: >> >> >> eqX = TransientTerm(var=phi) == ExplicitDiffusionTerm(coeff=D, var=phi) >> for step in range(steps): >> eqX.solve(dt=timeStepDuration) >> >> >> there would be alot of different terms. Where can I find an explanation >> on how to change these variables you mentioned into an equation one >> requires to run on fipy? >> >> Sergio >> >> >> Sergio Manzetti >> >> <http://www.fjordforsk.no/logo_hr2.jpg> >> >> Fjordforsk AS <http://www.fjordforsk.no> >> <http://www.fjordforsk.no> >> Midtun >> 6894 Vangsnes >> Norge >> Org.nr. 911 659 654 >> Tlf: +47 57695621 <+47%2057%2069%2056%2021> >> Økolab <http://www.oekolab.com> | Nanofactory <http://www.nanofact.no> >> | AQ-Lab <http://www.aq-lab.no> | FAP <http://www.phap.no> >> >> >> ------------------------------ >> *From: *"Thibault Bridel-Bertomeu" <thibault.bridellel...@gmail.com> >> *To: *"fipy" <fipy@nist.gov> >> *Sent: *Thursday, May 18, 2017 2:09:47 PM >> *Subject: *Re: Problems running a simple PDE >> >> Hello again Sergio, >> When you define eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D) >> you ready the equation, and then eqX.solve(var=phi, dt=timeStepDuration) >> solves the equation for an infinitesimal timestep dt, and the variable >> phi as asked. >> >> You could also have written :: >> >> eqX = TransientTerm(var=phi) == ExplicitDiffusionTerm(coeff=D, var=phi) >> for step in range(steps): >> eqX.solve(dt=timeStepDuration) >> >> Best >> >> Thibault Bridel-Bertomeu >> >> >> 2017-05-18 13:57 GMT+02:00 Sergio Manzetti <sergio.manze...@fjordforsk.no >> >: >> >>> Dear all, referring to the email below, I removed the phi_nw and >>> phi_old, however, I am not sure where the actual PDE given in the first >>> example at . 58, >>> >>> d phi/dt = A d^2x∕dx^2 >>> >>> >>> appears in the actual python code. >>> >>> >>> This would be important to change the PDE into other PDEs >>> >>> Where do I define this PDE? >>> >>> Thanks >>> >>> Sergio Manzetti >>> >>> <http://www.fjordforsk.no/logo_hr2.jpg> >>> >>> Fjordforsk AS <http://www.fjordforsk.no> >>> <http://www.fjordforsk.no> >>> Midtun >>> 6894 Vangsnes >>> Norge >>> Org.nr. 911 659 654 >>> Tlf: +47 57695621 <+47%2057%2069%2056%2021> >>> Økolab <http://www.oekolab.com> | Nanofactory >>> <http://www.nanofact.no> | AQ-Lab <http://www.aq-lab.no> | FAP >>> <http://www.phap.no> >>> >>> >>> ------------------------------ >>> *From: *"sergio manzetti" <sergio.manze...@fjordforsk.no> >>> *To: *"fipy" <fipy@nist.gov> >>> *Sent: *Thursday, May 18, 2017 1:36:57 PM >>> *Subject: *Problems running a simple PDE >>> >>> Hello, following the manual, at the first example with the 1D diffusion, >>> I have tried to make the python script for this (python2.7) and I get a >>> missing name, which is not defined in the example at all. >>> >>> This is the name "cell", which appears on p 58 in the manual. >>> >>> Please see this code which describes this example: >>> >>> >>> Thanks! >>> >>> >>> >>> #Python script for solving a Partial Differential Equation in a >>> diffusion case >>> >>> >>> >>> from fipy import * >>> nx = 50 >>> dx = 1 >>> mesh = Grid1D(nx=nx, dx=dx) >>> phi = CellVariable(name="Solution variable", >>> mesh=mesh, >>> value=0.) >>> >>> >>> #We'll use the coefficient set to D=1 >>> >>> D=1 >>> >>> # We give then 2000 time-steps for the computation, which is defined by >>> 90 percent of the maximum stable timestep, which is given >>> >>> timeStepDuration = 0.9 * dx**2 / (2 * D) >>> steps = 2000 >>> >>> >>> #Then we define the boundary conditions >>> >>> valueLeft = 1 >>> >>> valueRight = 0 >>> >>> #The boundary conditions are represented as faces around the exterior of >>> the mesh, so we define the constraint by the given values on the left and >>> right side of the boundary using the phi.constrain() command >>> >>> phi.constrain(valueLeft, mesh.facesRight) >>> phi.constrain(valueRight, mesh.facesLeft) >>> >>> # We can also omit giving boundary conditions, then the default bc is >>> equivalent to a zero gradient. >>> >>> #At this stage we define the partial differential equation as a >>> numerical problem to be solved >>> >>> for step in range(steps): >>> for j in range(cells): >>> phi_new[j] = phi_old[j] \ >>> + (D * dt / dx**2) * (phi_old[j+1] - 2 * phi_old[j] + >>> phi_old[j-1]) >>> time += dt >>> >>> #and additional code for the boundary conditions >>> >>> eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D) >>> >>> >>> #We then want to view the results of the calculation and use the command >>> Viewer() for this >>> >>> phiAnalytical = CellVariable(name="analytical value", >>> mesh=mesh) >>> if __name__ == '__main__': >>> viewer = Viewer(vars=(phi, phiAnalytical), >>> datamin=0., datamax=1.) >>> viewer.plot() >>> >>> #If we have a semi-infinite domain, then the solution for this PDE is >>> phi=1-erf(x/2(Dt)^(0.5)). This requires to import SciPy library, so we >>> import that. >>> >>> x = mesh.cellCenters[0] >>> t = timeStepDuration * steps >>> >>> try: >>> from scipy.special import erf >>> phiAnalytical.setValue(1-erf(x / (2 * numerix.sqrt(D * t)))) >>> except ImportError: >>> print ("The Scipy Library is not avaliable to test the solution to >>> the given trasient diffusion equation") >>> >>> # The equation is then solved by repeatedly looping >>> >>> for step in range(steps): >>> eqX.solve(var=phi, >>> dt=timeStepDuration) >>> if __name__ == '__main__': >>> viewer.plot() >>> print (phi.allclose(phiAnalytical, atol = 7e-4)) >>> 1 >>> >>> if __name__ == '__main__': >>> raw_input("Explicit transient diffusion. Press <return> to >>> proceed...") >>> >>> >>> _______________________________________________ >>> 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 ]
