Re: [fipy] How do I get internal boundary conditions? letting one square remain within another.

2021-03-05 Thread VICTOR PEREZ 
Thanks, it worked for me.

El viernes, 5 de marzo de 2021 a las 11:19:08 UTC-5, [email protected] 
escribió:

> Victor - 
>
> An internal constraint, as explained at
>
>
> https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#internal-fixed-value
>
> is applied by a combination of implicit and explicit sources, not 
> diffusion terms, e.g.,
>
> eq = (TransientTerm() == DiffusionTerm(coeff=alpha0)
>   - ImplicitSourceTerm(coeff=mask*largeValue) 
>   + mask * largeValue * 30.)
>
>
> The idea is that, where `mask` is True, this equation effectively becomes
>
> eq = (0 == - ImplicitSourceTerm(coeff=mask*largeValue) + mask * largeValue 
> * 30.)
>
>
> or
>
> eq = (ImplicitSourceTerm(coeff=mask) == mask * 30.)
>
> which is an implicit way to write
>
>
> T = 30
>
>
> when `mask` is True. 
>
> When `mask` is False, then these sources go away and the transient 
> diffusion equation dominates.
>
>
> A few other changes are needed to get this working:
>
> - Although `mesh.cellCenters` is a rank-1 `CellVariable` and 
> `mesh.cellCenters[0]` is a rank-0 `CellVariable`, unfortunately, 
> [tuple-unpacking](
> https://www.w3schools.com/python/python_tuples_unpack.asp) does *not* 
> result in x, y, and z containing `CellVariable` objects. This confuses FiPY 
> later. Instead, you can write
>
> x, y, z = mesh.cellCenters[0], mesh.cellCenters[1], mesh.cellCenters[2]
>
> - You need a mask. To mask out the middle 1/3 cube of your domain, I would 
> write
>
>   mask = ((x > L1 / 3) & (x < 2 * L1 / 3)
>   & (y > L2 / 3) & (y < 2 * L2 / 3)
>   & (z > L3 / 3) & (z < 2 * L3 / 3))
>
> - Your time step is not consistent with your diffusivity, I would try
>
>   timeStepDuration = 10*0.9*dx**2/(2.*alpha0)
>
>
> Putting it all together:
>
> ```
> #!/usr/bin/env python3
> '''
> include lib
> '''
> from fipy import *
> import time
>
> nx = ny = 30.
> nz = 6.
> dx = dy = dz = 1.
> L1 = dx*nx
> L2 = dy*ny
> L3 = dz*nz
>
> mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)
> x, y, z = mesh.cellCenters[0], mesh.cellCenters[1], mesh.cellCenters[2]
>
> mask = ((x > L1 / 3) & (x < 2 * L1 / 3)
> & (y > L2 / 3) & (y < 2 * L2 / 3)
> & (z > L3 / 3) & (z < 2 * L3 / 3))
>
> T = CellVariable(mesh=mesh, name = r"TEMPERATURE", value=23.0)
> #T[mask] = 30.
> '''
> value_constranin = [24, 24]
> value_constranin_C1 = value_constranin[1]
> print(value_constranin_C1)
> '''
> value_C1, value_C2 = 23, 23
>
> X, Y, Z = mesh.faceCenters
>
> C1 = ((mesh.facesLeft & (Y < L1))
> | (mesh.facesTop & (X < L1)))
>
> C2 = ((mesh.facesRight & (Y < L2))
> | (mesh.facesBottom & (X < L2)))
>
> T.constrain(value_C1, C1)
> T.constrain(value_C2, C2)
> T.constrain(60, mesh.facesBack)
> '''
> START: Propiedades termicas del material (Ceramica)
> '''
> k = 1.3 # [W/m·K] [J/m.K.s] conductividad termica
> rho = 2300 # [kg/m**3] desidad
> cp = 840 # [J/(kg·K)] calor especifico
>
> alpha0 = k/rho/cp#[m2/s] difusividad termica
> print(alpha0)
> alpha = FaceVariable(mesh=mesh, value=alpha0)
>
> '''
> END: Propiedades termicas del material (Ceramica)
> '''
> largeValue = 1e+7
> print("Alpha: ",alpha*largeValue)
> D = 0.1
> eq = (TransientTerm() == DiffusionTerm(coeff=alpha)
>   - ImplicitSourceTerm(mask*largeValue)
>   + mask*largeValue*30.)
>
> '''
> START: Leyenda de la grafica
> '''
> viewer = MayaviClient(vars=T,
> limits={'T min': 20, ' T max': 70},
> datamin=23, datamax=70,
> title="TEMPERATURE")
> '''
> END: Leyenda de la grafica
> '''
> from builtins import range
> timeStepDuration = 10*0.9*dx**2/(2.*alpha0)
> steps = 10
> t = 0.
> start_time = time.time()
> #print("start: ", start_time)
> print("start: t1 = ", timeStepDuration)
> for step in range(steps):
> eq.solve(var=T, dt=timeStepDuration)
> t = timeStepDuration + t
> #print("Time: ",step, time.time() - start_time)
> print("Time: t ", t)
> viewer.plot()
> ```
>
> On Mar 4, 2021, at 8:02 PM, VICTOR PEREZ  wrote:
>
> for the moment I have the conditions of the edge of the square. but the 
> intense conditions could not. In the image I show how the internal border 
> conditions should be.
>
> #!/usr/bin/env python3
> '''
> include lib
> '''
> from fipy import *
> import time
>
> nx = ny = 30.
> nz = 6.
> dx = dy = dz = 1.
> L1 = dx*nx
> L2 = dy*ny
> L3 = dz*nz
>
> mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)
> x, y, z = mesh.cellCenters
>
> #mask = mesh.cellCenters
> #mask = mesh.facesBack == 0.5
> #print("mask: ", mask)
> #mask0 = mesh.cellCenters.value[0] > True
> #mask1 = mesh.cellCenters.value[1] 
> #mask2 = mesh.cellCenters.value[2]
>
> #print("mask0: ", mask0)
> #print("mask1: ", mask1)
> #print("mask2: ", mask2)
>
> T = CellVariable(mesh=mesh, name = r"TEMPERATURE", value=23.0)
> #T[mask] = 30.
> '''
> value_constranin = [24, 24]
> value_constranin_C1 = value_constranin[1]
> print(value_constranin_C1)
> '''
> value_C1, value_C2 = 23, 23
>
> X, Y, Z = mesh.faceCenters
>
> C1 = ((mesh.facesLeft & (Y < L1)) 
> | (mesh.facesTop &

Re: [fipy] How do I get internal boundary conditions? letting one square remain within another.

2021-03-05 Thread 'Guyer, Jonathan E. Dr. (Fed)' via fipy 
Victor -

An internal constraint, as explained at

   https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#internal-fixed-value

is applied by a combination of implicit and explicit sources, not diffusion 
terms, e.g.,

eq = (TransientTerm() == DiffusionTerm(coeff=alpha0)
  - ImplicitSourceTerm(coeff=mask*largeValue)
  + mask * largeValue * 30.)

The idea is that, where `mask` is True, this equation effectively becomes

eq = (0 == - ImplicitSourceTerm(coeff=mask*largeValue) + mask * largeValue * 
30.)

or

eq = (ImplicitSourceTerm(coeff=mask) == mask * 30.)

which is an implicit way to write

T = 30

when `mask` is True.

When `mask` is False, then these sources go away and the transient diffusion 
equation dominates.


A few other changes are needed to get this working:

- Although `mesh.cellCenters` is a rank-1 `CellVariable` and 
`mesh.cellCenters[0]` is a rank-0 `CellVariable`, unfortunately, 
[tuple-unpacking](https://www.w3schools.com/python/python_tuples_unpack.asp) 
does *not* result in x, y, and z containing `CellVariable` objects. This 
confuses FiPY later. Instead, you can write

x, y, z = mesh.cellCenters[0], mesh.cellCenters[1], mesh.cellCenters[2]

- You need a mask. To mask out the middle 1/3 cube of your domain, I would write

  mask = ((x > L1 / 3) & (x < 2 * L1 / 3)
  & (y > L2 / 3) & (y < 2 * L2 / 3)
  & (z > L3 / 3) & (z < 2 * L3 / 3))

- Your time step is not consistent with your diffusivity, I would try

  timeStepDuration = 10*0.9*dx**2/(2.*alpha0)


Putting it all together:

```
#!/usr/bin/env python3
'''
include lib
'''
from fipy import *
import time

nx = ny = 30.
nz = 6.
dx = dy = dz = 1.
L1 = dx*nx
L2 = dy*ny
L3 = dz*nz

mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)
x, y, z = mesh.cellCenters[0], mesh.cellCenters[1], mesh.cellCenters[2]

mask = ((x > L1 / 3) & (x < 2 * L1 / 3)
& (y > L2 / 3) & (y < 2 * L2 / 3)
& (z > L3 / 3) & (z < 2 * L3 / 3))

T = CellVariable(mesh=mesh, name = r"TEMPERATURE", value=23.0)
#T[mask] = 30.
'''
value_constranin = [24, 24]
value_constranin_C1 = value_constranin[1]
print(value_constranin_C1)
'''
value_C1, value_C2 = 23, 23

X, Y, Z = mesh.faceCenters

C1 = ((mesh.facesLeft & (Y < L1))
| (mesh.facesTop & (X < L1)))

C2 = ((mesh.facesRight & (Y < L2))
| (mesh.facesBottom & (X < L2)))

T.constrain(value_C1, C1)
T.constrain(value_C2, C2)
T.constrain(60, mesh.facesBack)
'''
START: Propiedades termicas del material (Ceramica)
'''
k = 1.3 # [W/m·K] [J/m.K.s] conductividad termica
rho = 2300 # [kg/m**3] desidad
cp = 840 # [J/(kg·K)] calor especifico

alpha0 = k/rho/cp#[m2/s] difusividad termica
print(alpha0)
alpha = FaceVariable(mesh=mesh, value=alpha0)

'''
END: Propiedades termicas del material (Ceramica)
'''
largeValue = 1e+7
print("Alpha: ",alpha*largeValue)
D = 0.1
eq = (TransientTerm() == DiffusionTerm(coeff=alpha)
  - ImplicitSourceTerm(mask*largeValue)
  + mask*largeValue*30.)

'''
START: Leyenda de la grafica
'''
viewer = MayaviClient(vars=T,
limits={'T min': 20, ' T max': 70},
datamin=23, datamax=70,
title="TEMPERATURE")
'''
END: Leyenda de la grafica
'''
from builtins import range
timeStepDuration = 10*0.9*dx**2/(2.*alpha0)
steps = 10
t = 0.
start_time = time.time()
#print("start: ", start_time)
print("start: t1 = ", timeStepDuration)
for step in range(steps):
eq.solve(var=T, dt=timeStepDuration)
t = timeStepDuration + t
#print("Time: ",step, time.time() - start_time)
print("Time: t ", t)
viewer.plot()
```

On Mar 4, 2021, at 8:02 PM, VICTOR PEREZ 
mailto:[email protected]>> wrote:

for the moment I have the conditions of the edge of the square. but the intense 
conditions could not. In the image I show how the internal border conditions 
should be.

#!/usr/bin/env python3
'''
include lib
'''
from fipy import *
import time

nx = ny = 30.
nz = 6.
dx = dy = dz = 1.
L1 = dx*nx
L2 = dy*ny
L3 = dz*nz

mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)
x, y, z = mesh.cellCenters

#mask = mesh.cellCenters
#mask = mesh.facesBack == 0.5
#print("mask: ", mask)
#mask0 = mesh.cellCenters.value[0] > True
#mask1 = mesh.cellCenters.value[1]
#mask2 = mesh.cellCenters.value[2]

#print("mask0: ", mask0)
#print("mask1: ", mask1)
#print("mask2: ", mask2)

T = CellVariable(mesh=mesh, name = r"TEMPERATURE", value=23.0)
#T[mask] = 30.
'''
value_constranin = [24, 24]
value_constranin_C1 = value_constranin[1]
print(value_constranin_C1)
'''
value_C1, value_C2 = 23, 23

X, Y, Z = mesh.faceCenters

C1 = ((mesh.facesLeft & (Y < L1))
| (mesh.facesTop & (X < L1)))

C2 = ((mesh.facesRight & (Y < L2))
| (mesh.facesBottom & (X < L2)))

T.constrain(value_C1, C1)
T.constrain(value_C2, C2)
T.constrain(60, mesh.facesBack)
'''
START: Propiedades termicas del material (Ceramica)
'''
k = 1.3 # [W/m·K] [J/m.K.s] conductividad termica
rho = 2300 # [kg/m**3] desidad
cp = 840 # [J/(kg·K)] calor especifico

alpha0 = k/rho/cp#[m2/s] difusividad termica
print(alpha0)
alpha = FaceVariable(mesh=mesh, value