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 <[email protected]>: > 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" <[email protected]> > *To: *"fipy" <[email protected]> > *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 <[email protected]> > : > >> 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" <[email protected]> >> *To: *"fipy" <[email protected]> >> *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 >> [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 ]
