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 ] [ http://www.fjordforsk.no/ | Fjordforsk AS ] [ http://www.fjordforsk.no/ | ] Midtun 6894 Vangsnes Norge Org.nr. 911 659 654 Tlf: +47 57695621 [ http://www.oekolab.com/ | Økolab ] | [ http://www.nanofact.no/ | Nanofactory ] | [ http://www.aq-lab.no/ | AQ-Lab ] | [ http://www.phap.no/ | FAP ] 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 | 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 | 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 | 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 < [ mailto:sergio.manze...@fjordforsk.no | 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 ] [ http://www.fjordforsk.no/ | Fjordforsk AS ] [ http://www.fjordforsk.no/ | ] Midtun 6894 Vangsnes Norge Org.nr. 911 659 654 Tlf: [ tel:+47%2057%2069%2056%2021 | +47 57695621 ] [ http://www.oekolab.com/ | Økolab ] | [ http://www.nanofact.no/ | Nanofactory ] | [ http://www.aq-lab.no/ | AQ-Lab ] | [ http://www.phap.no/ | FAP ] From: "Thibault Bridel-Bertomeu" < [ mailto:thibault.bridellel...@gmail.com | thibault.bridellel...@gmail.com ] > To: "fipy" < [ mailto:fipy@nist.gov | 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 < [ mailto:sergio.manze...@fjordforsk.no | sergio.manze...@fjordforsk.no ] > : BQ_BEGIN 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 ] [ http://www.fjordforsk.no/ | Fjordforsk AS ] [ http://www.fjordforsk.no/ | ] Midtun 6894 Vangsnes Norge Org.nr. 911 659 654 Tlf: [ tel:+47%2057%2069%2056%2021 | +47 57695621 ] [ http://www.oekolab.com/ | Økolab ] | [ http://www.nanofact.no/ | Nanofactory ] | [ http://www.aq-lab.no/ | AQ-Lab ] | [ http://www.phap.no/ | FAP ] From: "sergio manzetti" < [ mailto:sergio.manze...@fjordforsk.no | sergio.manze...@fjordforsk.no ] > To: "fipy" < [ mailto:fipy@nist.gov | 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 [ mailto:fipy@nist.gov | fipy@nist.gov ] [ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ] [ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy | https://email.nist.gov/mailman/listinfo/fipy ] ] _______________________________________________ fipy mailing list [ mailto:fipy@nist.gov | fipy@nist.gov ] [ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ] [ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy | https://email.nist.gov/mailman/listinfo/fipy ] ] _______________________________________________ fipy mailing list [ mailto:fipy@nist.gov | fipy@nist.gov ] [ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ] [ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy | https://email.nist.gov/mailman/listinfo/fipy ] ] BQ_END _______________________________________________ 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 ]
