Re: 1D diffusion problem
Hi Jon, Thank you so much for taking the time to explain everything! Clara > On 22 Apr 2019, at 13:31, Guyer, Jonathan E. Dr. (Fed) via fipy > wrote: > > Also, the natural boundary condition for FiPy is no-flux, so there's no need > to do > > T_var.faceGrad.constrain(0, where=mesh.facesLeft) > >> On Apr 22, 2019, at 1:28 PM, Guyer, Jonathan E. Dr. (Fed) via fipy >> wrote: >> >> Clara - >> >> - The cause of everything dropping to zero is that T_var.old is initialized >> to zero. This isn't a particularly useful initial value, but that's what it >> is. Usually, T_var.updateOld() is the first thing called in the time step >> loop, so this initial condition doesn't matter, but because of the way >> you've structured your adaptive time steps, the first sweep is done with an >> erroneous initial condition. >> >> Call T_var.updateOld() one time before you start taking time steps. >> >> - Your sweep tolerance is too severe. I adjusted to 1.0e-6 and the adaptive >> stepper works much better. >> >> - The way you define and update dt, you'll get recursion errors eventually. >> Instead, define >> >> dt = 0.9*dr**2/(2.0*Diff[0].value) >> >> If you're curious, this is because: >> >> * Diff is a FaceVariable >> * so, Diff[0] is a UnaryOperator Variable >> * so, dt = 0.9*dr**2/(2.0*Diff[0]) is a BinaryOperator Variable >> * so, dt *= 1.1/dt *= 0.8 is a BinaryOperator Variable that gets nested >> deeper >>and deeper every time step >> >> Diff[0].value short circuits that >> >> - Your problem is completely linear, so sweeps aren't necessary at this point >> >> - This equation is fully implicit, so there's no need to start with such a >> tiny initial time step >> >> - Jon >> >>> On Apr 19, 2019, at 2:42 PM, Clara Maurel wrote: >>> >>> Hello everyone, >>> >>> I apologize in advance for this very basic question… >>> I am trying to solve a 1D diffusion problem where I have a cooling planet >>> (core+mantle) with an initial temperature profile and a diffusion >>> coefficient that takes different values in the core and in the mantle. I >>> want a fixed boundary condition at the surface (T = 200 K) and a zero >>> gradient at the center of the body, but the temperature at the center (left >>> of the mesh) is not fixed and should cool. >>> When I run the script, the whole profile immediately goes down to zero >>> (except the end right of the mesh) and diffuses from there. I must be doing >>> something wrong with the left boundary condition but I can’t figure out >>> what. >>> >>> Could someone help me to find out what I missed? Any hint would be greatly >>> appreciated! >>> >>> Thank you in advance for your help! >>> Clara >>> >>> Here is my simple script: >>> >>> from datetime import datetime >>> startTime = datetime.now() >>> import numpy as np >>> import scipy >>> from fipy import * >>> from fipy.solvers.pysparse import LinearLUSolver as Solver >>> from sympy import * >>> from scipy import interpolate >>> >>> def Get_closest_id(L,value): >>> return list(L).index(min(L, key=lambda x:abs(x-value))) >>> >>> ## >>> # DIFFUSION EQUATION: INITIALIZATION # >>> ## >>> >>> ## Dimensions (km) >>> R_body = 200.0 >>> R_core = 50.0 >>> >>> ## Temperatures (K) >>> T_core = 1200.0 >>> T_surf = 200.0 >>> >>> ## Mesh >>> nr = 1000 >>> dr = R_body / nr >>> mesh = Grid1D(nx=nr, dx=dr) >>> >>> ## Variable >>> T_var = CellVariable(name="T", mesh=mesh, hasOld=True) >>> >>> ## Initial temperature profile >>> slope = (T_core-T_surf)/(R_core-R_body) >>> intercept = T_core-slope*R_core >>> T_var[:] = np.array([T_core]*int(R_core/dr) + >>> list(slope*np.arange(R_core+dr, R_body+dr, dr)+intercept)) >>> >>> ## Initial conditions (fixed surface value, null gradient at the center of >>> the core) >>> T_var.constrain(T_surf, mesh.facesRight) >>> T_var.faceGrad.constrain(0, where=mesh.facesLeft) >>> >>> ## Diffusion coefficient (different in core and mantle, in km2 My-1) >>> Diff = FaceVariable(mesh=mesh) >>> X = mesh.faceCenter
Re: 1D diffusion problem
Also, the natural boundary condition for FiPy is no-flux, so there's no need to
do
T_var.faceGrad.constrain(0, where=mesh.facesLeft)
> On Apr 22, 2019, at 1:28 PM, Guyer, Jonathan E. Dr. (Fed) via fipy
> wrote:
>
> Clara -
>
> - The cause of everything dropping to zero is that T_var.old is initialized
> to zero. This isn't a particularly useful initial value, but that's what it
> is. Usually, T_var.updateOld() is the first thing called in the time step
> loop, so this initial condition doesn't matter, but because of the way you've
> structured your adaptive time steps, the first sweep is done with an
> erroneous initial condition.
>
> Call T_var.updateOld() one time before you start taking time steps.
>
> - Your sweep tolerance is too severe. I adjusted to 1.0e-6 and the adaptive
> stepper works much better.
>
> - The way you define and update dt, you'll get recursion errors eventually.
> Instead, define
>
> dt = 0.9*dr**2/(2.0*Diff[0].value)
>
> If you're curious, this is because:
>
> * Diff is a FaceVariable
> * so, Diff[0] is a UnaryOperator Variable
> * so, dt = 0.9*dr**2/(2.0*Diff[0]) is a BinaryOperator Variable
> * so, dt *= 1.1/dt *= 0.8 is a BinaryOperator Variable that gets nested
> deeper
> and deeper every time step
>
> Diff[0].value short circuits that
>
> - Your problem is completely linear, so sweeps aren't necessary at this point
>
> - This equation is fully implicit, so there's no need to start with such a
> tiny initial time step
>
> - Jon
>
>> On Apr 19, 2019, at 2:42 PM, Clara Maurel wrote:
>>
>> Hello everyone,
>>
>> I apologize in advance for this very basic question…
>> I am trying to solve a 1D diffusion problem where I have a cooling planet
>> (core+mantle) with an initial temperature profile and a diffusion
>> coefficient that takes different values in the core and in the mantle. I
>> want a fixed boundary condition at the surface (T = 200 K) and a zero
>> gradient at the center of the body, but the temperature at the center (left
>> of the mesh) is not fixed and should cool.
>> When I run the script, the whole profile immediately goes down to zero
>> (except the end right of the mesh) and diffuses from there. I must be doing
>> something wrong with the left boundary condition but I can’t figure out what.
>>
>> Could someone help me to find out what I missed? Any hint would be greatly
>> appreciated!
>>
>> Thank you in advance for your help!
>> Clara
>>
>> Here is my simple script:
>>
>> from datetime import datetime
>> startTime = datetime.now()
>> import numpy as np
>> import scipy
>> from fipy import *
>> from fipy.solvers.pysparse import LinearLUSolver as Solver
>> from sympy import *
>> from scipy import interpolate
>>
>> def Get_closest_id(L,value):
>>return list(L).index(min(L, key=lambda x:abs(x-value)))
>>
>> ##
>> # DIFFUSION EQUATION: INITIALIZATION #
>> ##
>>
>> ## Dimensions (km)
>> R_body = 200.0
>> R_core = 50.0
>>
>> ## Temperatures (K)
>> T_core = 1200.0
>> T_surf = 200.0
>>
>> ## Mesh
>> nr = 1000
>> dr = R_body / nr
>> mesh = Grid1D(nx=nr, dx=dr)
>>
>> ## Variable
>> T_var = CellVariable(name="T", mesh=mesh, hasOld=True)
>>
>> ## Initial temperature profile
>> slope = (T_core-T_surf)/(R_core-R_body)
>> intercept = T_core-slope*R_core
>> T_var[:] = np.array([T_core]*int(R_core/dr) +
>> list(slope*np.arange(R_core+dr, R_body+dr, dr)+intercept))
>>
>> ## Initial conditions (fixed surface value, null gradient at the center of
>> the core)
>> T_var.constrain(T_surf, mesh.facesRight)
>> T_var.faceGrad.constrain(0, where=mesh.facesLeft)
>>
>> ## Diffusion coefficient (different in core and mantle, in km2 My-1)
>> Diff = FaceVariable(mesh=mesh)
>> X = mesh.faceCenters[0]
>> Diff.setValue(150.0, where=(R_core >= X))
>> Diff.setValue(15.0, where=(R_core < X) & (R_body >= X))
>>
>> ## Equation
>> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff)
>>
>> if __name__ == "__main__":
>>viewer = Viewer(vars=(T_var),ymin=0,ymax=1400)
>>viewer.plot()
>>
>> if __name__ == '__main__':
>>raw_input("Press to proceed...")
>>
>> ## Time and integration parameters (time in My; 3.154e13 = 1 My in s)
>> time = 0.0
>> duration = 1000.0
>> dt = 0.9*dr**2/(2.0
Re: 1D diffusion problem
Clara -
- The cause of everything dropping to zero is that T_var.old is initialized to
zero. This isn't a particularly useful initial value, but that's what it is.
Usually, T_var.updateOld() is the first thing called in the time step loop, so
this initial condition doesn't matter, but because of the way you've structured
your adaptive time steps, the first sweep is done with an erroneous initial
condition.
Call T_var.updateOld() one time before you start taking time steps.
- Your sweep tolerance is too severe. I adjusted to 1.0e-6 and the adaptive
stepper works much better.
- The way you define and update dt, you'll get recursion errors eventually.
Instead, define
dt = 0.9*dr**2/(2.0*Diff[0].value)
If you're curious, this is because:
* Diff is a FaceVariable
* so, Diff[0] is a UnaryOperator Variable
* so, dt = 0.9*dr**2/(2.0*Diff[0]) is a BinaryOperator Variable
* so, dt *= 1.1/dt *= 0.8 is a BinaryOperator Variable that gets nested
deeper
and deeper every time step
Diff[0].value short circuits that
- Your problem is completely linear, so sweeps aren't necessary at this point
- This equation is fully implicit, so there's no need to start with such a tiny
initial time step
- Jon
> On Apr 19, 2019, at 2:42 PM, Clara Maurel wrote:
>
> Hello everyone,
>
> I apologize in advance for this very basic question…
> I am trying to solve a 1D diffusion problem where I have a cooling planet
> (core+mantle) with an initial temperature profile and a diffusion coefficient
> that takes different values in the core and in the mantle. I want a fixed
> boundary condition at the surface (T = 200 K) and a zero gradient at the
> center of the body, but the temperature at the center (left of the mesh) is
> not fixed and should cool.
> When I run the script, the whole profile immediately goes down to zero
> (except the end right of the mesh) and diffuses from there. I must be doing
> something wrong with the left boundary condition but I can’t figure out what.
>
> Could someone help me to find out what I missed? Any hint would be greatly
> appreciated!
>
> Thank you in advance for your help!
> Clara
>
> Here is my simple script:
>
> from datetime import datetime
> startTime = datetime.now()
> import numpy as np
> import scipy
> from fipy import *
> from fipy.solvers.pysparse import LinearLUSolver as Solver
> from sympy import *
> from scipy import interpolate
>
> def Get_closest_id(L,value):
> return list(L).index(min(L, key=lambda x:abs(x-value)))
>
> ##
> # DIFFUSION EQUATION: INITIALIZATION #
> ##
>
> ## Dimensions (km)
> R_body = 200.0
> R_core = 50.0
>
> ## Temperatures (K)
> T_core = 1200.0
> T_surf = 200.0
>
> ## Mesh
> nr = 1000
> dr = R_body / nr
> mesh = Grid1D(nx=nr, dx=dr)
>
> ## Variable
> T_var = CellVariable(name="T", mesh=mesh, hasOld=True)
>
> ## Initial temperature profile
> slope = (T_core-T_surf)/(R_core-R_body)
> intercept = T_core-slope*R_core
> T_var[:] = np.array([T_core]*int(R_core/dr) + list(slope*np.arange(R_core+dr,
> R_body+dr, dr)+intercept))
>
> ## Initial conditions (fixed surface value, null gradient at the center of
> the core)
> T_var.constrain(T_surf, mesh.facesRight)
> T_var.faceGrad.constrain(0, where=mesh.facesLeft)
>
> ## Diffusion coefficient (different in core and mantle, in km2 My-1)
> Diff = FaceVariable(mesh=mesh)
> X = mesh.faceCenters[0]
> Diff.setValue(150.0, where=(R_core >= X))
> Diff.setValue(15.0, where=(R_core < X) & (R_body >= X))
>
> ## Equation
> eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff)
>
> if __name__ == "__main__":
> viewer = Viewer(vars=(T_var),ymin=0,ymax=1400)
> viewer.plot()
>
> if __name__ == '__main__':
> raw_input("Press to proceed...")
>
> ## Time and integration parameters (time in My; 3.154e13 = 1 My in s)
> time = 0.0
> duration = 1000.0
> dt = 0.9*dr**2/(2.0*Diff[0])
>
> total_sweeps = 4
> tolerance = 1.0e-10
>
> solver = Solver()
>
> #---#
> # DIFFUSION EQUATION: MAIN LOOP #
> #---#
>
> while time < duration:
>
> res0 = eq.sweep(T_var, dt=dt, solver=solver)
>
> for sweeps in range(total_sweeps):
> res = eq.sweep(T_var, dt=dt, solver=solver)
>
> if res < res0 * tolerance:
> time += dt
> dt *= 1.1
> T_var.updateOld()
> print dt, res
> if __name__ == '__main__':
> viewer.plot()
>
> else:
> dt *= 0.8
> T_var[:] = T_var.old
> print "w
1D diffusion problem
Hello everyone,
I apologize in advance for this very basic question…
I am trying to solve a 1D diffusion problem where I have a cooling planet
(core+mantle) with an initial temperature profile and a diffusion coefficient
that takes different values in the core and in the mantle. I want a fixed
boundary condition at the surface (T = 200 K) and a zero gradient at the center
of the body, but the temperature at the center (left of the mesh) is not fixed
and should cool.
When I run the script, the whole profile immediately goes down to zero (except
the end right of the mesh) and diffuses from there. I must be doing something
wrong with the left boundary condition but I can’t figure out what.
Could someone help me to find out what I missed? Any hint would be greatly
appreciated!
Thank you in advance for your help!
Clara
Here is my simple script:
from datetime import datetime
startTime = datetime.now()
import numpy as np
import scipy
from fipy import *
from fipy.solvers.pysparse import LinearLUSolver as Solver
from sympy import *
from scipy import interpolate
def Get_closest_id(L,value):
return list(L).index(min(L, key=lambda x:abs(x-value)))
##
# DIFFUSION EQUATION: INITIALIZATION #
##
## Dimensions (km)
R_body = 200.0
R_core = 50.0
## Temperatures (K)
T_core = 1200.0
T_surf = 200.0
## Mesh
nr = 1000
dr = R_body / nr
mesh = Grid1D(nx=nr, dx=dr)
## Variable
T_var = CellVariable(name="T", mesh=mesh, hasOld=True)
## Initial temperature profile
slope = (T_core-T_surf)/(R_core-R_body)
intercept = T_core-slope*R_core
T_var[:] = np.array([T_core]*int(R_core/dr) + list(slope*np.arange(R_core+dr,
R_body+dr, dr)+intercept))
## Initial conditions (fixed surface value, null gradient at the center of the
core)
T_var.constrain(T_surf, mesh.facesRight)
T_var.faceGrad.constrain(0, where=mesh.facesLeft)
## Diffusion coefficient (different in core and mantle, in km2 My-1)
Diff = FaceVariable(mesh=mesh)
X = mesh.faceCenters[0]
Diff.setValue(150.0, where=(R_core >= X))
Diff.setValue(15.0, where=(R_core < X) & (R_body >= X))
## Equation
eq = TransientTerm(coeff=1.) == DiffusionTerm(Diff)
if __name__ == "__main__":
viewer = Viewer(vars=(T_var),ymin=0,ymax=1400)
viewer.plot()
if __name__ == '__main__':
raw_input("Press to proceed...")
## Time and integration parameters (time in My; 3.154e13 = 1 My in s)
time = 0.0
duration = 1000.0
dt = 0.9*dr**2/(2.0*Diff[0])
total_sweeps = 4
tolerance = 1.0e-10
solver = Solver()
#---#
# DIFFUSION EQUATION: MAIN LOOP #
#---#
while time < duration:
res0 = eq.sweep(T_var, dt=dt, solver=solver)
for sweeps in range(total_sweeps):
res = eq.sweep(T_var, dt=dt, solver=solver)
if res < res0 * tolerance:
time += dt
dt *= 1.1
T_var.updateOld()
print dt, res
if __name__ == '__main__':
viewer.plot()
else:
dt *= 0.8
T_var[:] = T_var.old
print "warning " + str(res)
print datetime.now() - startTime
if __name__ == '__main__':
raw_input("Press to proceed...")
smime.p7s
Description: S/MIME cryptographic signature
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Using FiPy to solve a 1D diffusion problem under a potential profile
Other than a change of sign on the ConvectionTerm, I don't see any problem with
how you wrote your equation.
There's not much to do to debug a seg fault, other than isolate what line of
code is causing it. Comment out all the lines in your code and successively add
them until the seg fault occurs.
How big is your mesh?
If you create a Viewer with U_z, does it look like you think it should from
your experimental data?
What are the four failures when you run `python setup.py test`? If something
deeply was wrong, I'd expect many more than 4 failures, but it still might be
helpful to know.
- Jon
> On Jul 18, 2017, at 4:42 AM, WANG Lingjian <wan...@dc.tohoku.ac.jp> wrote:
>
> Hello everyone, I am trying to solve a 1D diffusion problem in a shape like
> the following equation using FiPy. Where U(z) is a 1D potential fitted from
> experiment data. I created this U_z as a CellVariable like the following: ```
> U_z = CellVariable(name="potential", mesh=mesh)
> U_z.setValue(fitted_function(x) * D * beta) ``` >From my understanding, this
> equation can be interpreted using FiPy syntax as ``` F_z = U_z.faceGrad eqX =
> TransientTerm() == (DiffusionTerm(coeff=D) + ConvectionTerm(coeff=F_z)) ```
> But when I try to solve it using commands that look like this, ``` for step
> in range(steps): eqX.solve(var=phi, dt=timeStepDuration) print "time step at
> {}".format(step) viewer.plot() ``` I encounter a Segmentation fault error
> with no further information indicating any possible error. Could anyone tell
> me: #1 Is there any misunderstanding in the above processes(i.e. treating the
> dU/dz term as a ConvectionTerm, using .faceGrad ...)? #2 Is there any option
> for debugging the !
Segmentation fault error? Thanks in advance. In case that software environment
issue, I provide some system information and FiPy test information below.
[SYSTEM INFO] System: fedora17(Beefy Miracle) Python version: anaconda2 4.4.0
FiPy installation method: "conda create --name MYFIPYENV --channel guyer
--channel conda-forge fipy" . [FiPy test status] Testing of FiPy using ```
python -c "import fipy; fipy.test()" ``` Fails with AttributeError: 'module'
object has no attribute 'testFiPy'. But test using ``` python setup.py test ```
under the FiPy source folder gives ### FAILED (failures=4)
!!!
Skipped 45 doctest examples because neither `lsmlib` nor `skfmm` can be found
on the $PATH Skipped 2 doctest examples because `lsmlib` must be used to run
some tests Skipped 1 doctest examples because `skfmm` must be used to run some
tests ### _!
__
> 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 ]
Using FiPy to solve a 1D diffusion problem under a potential profile
Hello everyone,
I am trying to solve a 1D diffusion problem in a shape like the following
equation using FiPy.
http://imgur.com/35xl6Cg;>http://i.imgur.com/35xl6Cg.png;
title="source: imgur.com" />
Where U(z) is a 1D potential fitted from experiment data. I created this U_z as
a CellVariable like the following:
```
U_z = CellVariable(name="potential", mesh=mesh)
U_z.setValue(fitted_function(x) * D * beta)
```
>From my understanding, this equation can be interpreted using FiPy syntax as
```
F_z = U_z.faceGrad
eqX = TransientTerm() == (DiffusionTerm(coeff=D) + ConvectionTerm(coeff=F_z))
```
But when I try to solve it using commands that look like this,
```
for step in range(steps):
eqX.solve(var=phi, dt=timeStepDuration)
print "time step at {}".format(step)
viewer.plot()
```
I encounter a Segmentation fault error with no further information indicating
any possible error.
Could anyone tell me:
#1 Is there any misunderstanding in the above processes(i.e. treating the dU/dz
term as a ConvectionTerm, using .faceGrad ...)?
#2 Is there any option for debugging the Segmentation fault error?
Thanks in advance.
In case that software environment issue, I provide some system information and
FiPy test information below.
[SYSTEM INFO]
System: fedora17(Beefy Miracle)
Python version: anaconda2 4.4.0
FiPy installation method: "conda create --name MYFIPYENV --channel guyer
--channel conda-forge fipy" .
[FiPy test status]
Testing of FiPy using
```
python -c "import fipy; fipy.test()"
```
Fails with AttributeError: 'module' object has no attribute 'testFiPy'.
But test using
```
python setup.py test
```
under the FiPy source folder gives
###
FAILED (failures=4)
!!!
Skipped 45 doctest examples because neither `lsmlib` nor `skfmm` can be found
on the $PATH
Skipped 2 doctest examples because `lsmlib` must be used to run some tests
Skipped 1 doctest examples because `skfmm` must be used to run some tests
###
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Flux=0 constrains on both sizes cause problems with simple 1D diffusion problem
Dear list, I have an issue with a simple script, (here https://gist.github.com/4037850 ). When I choose flux=0 boundary conditions on both size of this 1D diffusion problem, I don't think the results are correct. Physically this describes something like the absorption and emission of light. With the flux=0 constraints I would expect the value of phi (in this case electron density) to build up inside the volume because it cannot escape. However, the fipy solution for this case doesn't look physical. I'm probably missing something! Best regards, Dan. PS. Latex version of the equations http://imgur.com/mCcX2 . ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Flux=0 constrains on both sizes cause problems with simple 1D diffusion problem
Hi Daniel, In your script, you are solving an steady state problem. Add a transient term, and preferably sweep to get the solution. Perhaps sth like this: equation = (TransientTerm() == diffusion_term + source_term) timestep = your time step steps = 10 for step in range(steps) res = 1.e6 phi.updateOld() while res1.e-5 res = equation.sweep(var = phi, dt = timestep) Kind regards, Ali On Thu, 2012-11-08 at 19:02 +0900, Daniel Farrell wrote: Dear list, I have an issue with a simple script, (here https://gist.github.com/4037850). When I choose flux=0 boundary conditions on both size of this 1D diffusion problem, I don't think the results are correct. Physically this describes something like the absorption and emission of light. With the flux=0 constraints I would expect the value of phi (in this case electron density) to build up inside the volume because it cannot escape. However, the fipy solution for this case doesn't look physical. I'm probably missing something! Best regards, Dan. PS. Latex version of the equations http://imgur.com/mCcX2 . ___ 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 ]
