A file ? You can export the result to a file yes, although I did not demonstrate it in my script. This is not fipy though, it’s about using numpy to write an array of values into a file !
As for the comparison with an analytic solution, I invite you to take a look at the mesh1D diffusion example where they create a phiAnalytical. To plot it, you can do vi = Viewer(vars=(phi, phiAnalytical)) vi.plot() which should plot at the same time the solution of the iterative process from fipy and the analytical value. Cordialement, -- T. BRIDEL-BERTOMEU 2017-05-18 14:46 GMT+02:00 Sergio Manzetti <sergio.manze...@fjordforsk.no>: > OK, I see now. Where is this file you mention accessed? It does not come > out in the working directory. > > > 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:47:24 PM > > *Subject: *Re: Problems running a simple PDE > > The PDE you want to solve has a transient term, which is why the function > evolves in time. > However, as you can see it eventually converges to a given shape and stay > in that shape : this « stabilized shape » corresponds to the steady-state > result of your equation. > All you can see before are the transient organizations. > > As I said in my previous e-mail, you can access phi as a numpy array and > do with it as you want : this phi is the discrete representation of your > steady-state result. > > Cordialement, > > -- > T. BRIDEL-BERTOMEU > > > > 2017-05-18 14:40 GMT+02:00 Sergio Manzetti <sergio.manze...@fjordforsk.no> > : > >> I am just not sure on how to understand the actual beautiful result from >> the PDE. >> >> When initializing is pressed, it animates into a wave-like pattern. >> >> I thought the numerical solution to the PDE was an actual function, but >> this animated plot must be some form of time-dependent changing function? >> How can it be given as a function in addition to the graph? >> >> >> 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 ] > > _______________________________________________ > 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 ]
