Thank you very much its working! What I m trying to do is simulate results
from COMSOL with FiPy, however although I have a small time step 2E-04, I
cannot make a perfect fit. So, I suppose its the mesh... Is there anyway
to optimize it without trial and error as there are many faces?
COMSOL FIPY 0.01 1.3106662 1.318291217 0.02 1.638119 1.65445279 0.03
1.8321507 1.85201841 0.04 1.9702077 1.992648552 0.05 2.0775485 2.101983606
0.06 2.1653385 2.191437313 0.07 2.2394807 2.267119244 0.08 2.3039784
2.332693914 0.09 2.360845 2.39053534 0.1 2.4116504 2.442270102 0.11
2.4575982 2.489061483 0.12 2.4995654 2.531770171 0.13 2.5381932
2.571050769 0.14 2.5739691 2.607412613 0.15 2.6072903 2.641259674
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from fipy import *
import numpy as np
import time
def drange(start, end, step):
while start <= end:
yield start
start += step
def calcDT(cellSize, wire, wireI, silicone, solid1, solid1I, solid2,
solid2I):
t0 = time.time()
valuesX = valuesXLog = valuesY = np.array([])
gmsh = '''
cellSize = ''' + str(cellSize) + ''';
wireradius = ''' + str(wire['r']) + ''';
solidwidth = ''' + str(solid1['width']) + ''';
solid1thick = ''' + str(solid1['thick']) + ''';
solid1ithick = ''' + str(solid1I['thick']) + ''';
solid2thick = ''' + str(solid2['thick']) + ''';
solid2ithick = ''' + str(solid2I['thick']) + ''';
siliconethick = ''' + str(silicone['thick']) + ''';
siliconeithick = ''' + str(wireI['thick']) + ''';
Point(1) = {0, 0, 0, cellSize};
/* Wire */
Point(2) = {wireradius, 0, 0, cellSize};
Point(3) = {0, wireradius, 0, cellSize};
Point(4) = {0, -wireradius, 0, cellSize};
/* Silicone Interface */
Point(5) = {wireradius+siliconeithick, 0, 0, cellSize*5};
Point(6) = {0, wireradius+siliconeithick, 0, cellSize*5};
Point(7) = {0, -(wireradius+siliconeithick), 0, cellSize*5};
/* Silicone */
Point(8) = {0, siliconethick/2, 0, cellSize*80};
Point(9) = {0, -siliconethick/2, 0, cellSize*80};
Point(10) = {solidwidth, siliconethick/2, 0, cellSize*200};
Point(11) = {solidwidth, -siliconethick/2, 0, cellSize*200};
/* Solid 1 Interface */
Point(12) = {0, siliconethick/2+solid1ithick, 0, cellSize*200};
Point(13) = {solidwidth, siliconethick/2+solid1ithick, 0, cellSize*200};
/* Solid 1 */
Point(14) = {0, siliconethick/2+solid1ithick+solid1thick, 0, cellSize*500};
Point(15) = {solidwidth, siliconethick/2+solid1ithick+solid1thick, 0,
cellSize*2000};
/* Solid 2 Interface */
Point(16) = {0, -(siliconethick/2+solid2ithick), 0, cellSize*200};
Point(17) = {solidwidth, -(siliconethick/2+solid2ithick), 0, cellSize*200};
/* Solid 2 */
Point(18) = {0, -(siliconethick/2+solid2ithick+solid2thick), 0,
cellSize*500};
Point(19) = {solidwidth, -(siliconethick/2+solid2ithick+solid2thick), 0,
cellSize*2000};
Line(1) = {14, 15};
Line(2) = {15, 13};
Line(3) = {13, 12};
Line(4) = {12, 14};
Line(5) = {12, 8};
Line(6) = {8, 10};
Line(7) = {10, 13};
Line(8) = {10, 11};
Line(9) = {11, 9};
Line(10) = {9, 7};
Line(11) = {7, 4};
Line(12) = {4, 1};
Line(13) = {1, 3};
Line(14) = {3, 6};
Line(15) = {6, 8};
Line(16) = {9, 16};
Line(17) = {16, 17};
Line(18) = {17, 11};
Line(19) = {17, 19};
Line(20) = {19, 18};
Line(21) = {18, 16};
Circle(25) = {3, 1, 2};
Circle(26) = {2, 1, 4};
Circle(27) = {7, 1, 5};
Circle(28) = {5, 1, 6};
Line Loop(29) = {26, 12, 13, 25};
Plane Surface(30) = {29};/*"Wire"*/
Line Loop(31) = {27, 28, -14, 25, 26, -11};
Plane Surface(32) = {31};/*"Wire Interface"*/
Line Loop(33) = {10, 27, 28, 15, 6, 8, 9};
Plane Surface(34) = {33};/*"Silicone"*/
Line Loop(35) = {6, 7, 3, 5};
Plane Surface(36) = {35};/*"Interface 1"*/
Line Loop(37) = {4, 1, 2, 3};
Plane Surface(38) = {37};/*"Solid 1"*/
Line Loop(39) = {9, 16, 17, 18};
Plane Surface(40) = {39};/*"Interface 2"*/
Line Loop(41) = {19, 20, 21, 17};
Plane Surface(42) = {41};/*"Solid 2"*/
Physical Surface("Wire") = {30};
Physical Surface("WireI") = {32};
Physical Surface("Silicone") = {34};
Physical Surface("Solid1I") = {36};
Physical Surface("Solid1") = {38};
Physical Surface("Solid2I") = {40};
Physical Surface("Solid2") = {42};
'''
# print gmsh
mesh = GmshImporter2D(gmsh)
x, y = mesh.cellCenters[0], mesh.cellCenters[1]
dT = CellVariable(name=r'$\Delta T$', mesh=mesh, value=0.)
QCell = CellVariable(name=r'$Q$', mesh=mesh, value=0.)
tcCell = CellVariable(name=r'$Tc$', mesh=mesh, value=0.)
dencpCell = CellVariable(name=r'$DenCp$', mesh=mesh, value=0.)
x0 = y0 = 0
# Set Wire Interface Properties
meshWireI = mesh.physicalCells["WireI"]
tcCell.setValue(wireI['tc'], where=meshWireI)
dencpCell.setValue(wireI['den'] * wireI['cp'], where=meshWireI)
# Set Wire Properties
meshWire = mesh.physicalCells["Wire"]
QCell.setValue(wire['Q'], where=meshWire)
tcCell.setValue(wire['tc'], where=meshWire)
dencpCell.setValue(wire['den'] * wire['cp'], where=meshWire)
# Set Silicone Properties
meshSilicone = mesh.physicalCells["Silicone"]
tcCell.setValue(silicone['tc'], where=meshSilicone)
dencpCell.setValue(silicone['den'] * silicone['cp'], where=meshSilicone)
# Set Solid 1 Interface Properties
meshSolid1I = mesh.physicalCells["Solid1I"]
tcCell.setValue(solid1I['tc'], where=meshSolid1I)
dencpCell.setValue(solid1I['den'] * solid1I['cp'], where=meshSolid1I)
# Set Solid 1 Properties
meshSolid1 = mesh.physicalCells["Solid1"]
tcCell.setValue(solid1['tc'], where=meshSolid1)
dencpCell.setValue(solid1['den'] * solid1['cp'], where=meshSolid1)
# Set Solid 2 Interface Properties
meshSolid2I = mesh.physicalCells["Solid2I"]
tcCell.setValue(solid2I['tc'], where=meshSolid2I)
dencpCell.setValue(solid2I['den'] * solid2I['cp'], where=meshSolid2I)
# Set Solid 2 Properties
meshSolid2 = mesh.physicalCells["Solid2"]
tcCell.setValue(solid2['tc'], where=meshSolid2)
dencpCell.setValue(solid2['den'] * solid2['cp'], where=meshSolid2)
# Room temp
X, Y = mesh.faceCenters
mask = ((mesh.facesRight & (X == solid1['width'])) | (mesh.facesTop & (Y ==
silicone['thick'] / 2 + solid1I['thick'] + solid1['thick'])) |
(mesh.facesBottom & (Y == -silicone['thick'] / 2 - solid2I['thick'] -
solid2['thick'])))
dT.constrain(0, where=mask)
heatEq = (TransientTerm(coeff=dencpCell) == DiffusionTerm(coeff=tcCell) +
QCell)
timeStep = 0.0002
elapsedTime = 0.0
i = 1;
while(elapsedTime < 0.5):
heatEq.solve(var=dT, dt=timeStep)
dTVal = dT(((wire['r'] / numerix.sqrt(2),), (wire['r'] /
numerix.sqrt(2),)), order=1)[0]
elapsedTime += timeStep
if(__name__ == '__main__' and i % (0.01 / timeStep) == 0.0): # and
(elapsedTime * 1000) % (timeStep * 1000 * 10) == 0 # and
elapsedTime%0.01==0 and elapsedTime in times)
print '%s\t%s' % (elapsedTime, dTVal)
i = i + 1
return {'x':valuesX, 'xlog':valuesXLog, 'y':valuesY }
if __name__ == '__main__':
cellSize = 2e-06 # 5E-07
wire = {'r':1.2915E-05, 'Q':2.041289613757E+09, 'tc':57.5, 'den':16600.19,
'cp':150.5 }
wireI = {'thick':1E-05, 'tc':0.16, 'den':1000, 'cp':1600}
silicone = {'thick':9E-04, 'tc':0.1764, 'den':1000, 'cp':1800}
solid1 = {'thick':2E-02, 'width':0.025, 'tc':3.76, 'den':2596, 'cp':765}
solid1I = {'thick':1E-04, 'tc':3.76, 'den':2596, 'cp':765}
solid2 = {'thick':2E-02, 'width':0.025, 'tc':0.188, 'den':1289, 'cp':700}
solid2I = {'thick':1E-04, 'tc':2.8, 'den':1289, 'cp':700}
result = calcDT(cellSize, wire, wireI, silicone, solid1, solid1I, solid2,
solid2I)
On Mon, Jun 10, 2013 at 11:45 PM, Daniel Wheeler
<[email protected]>wrote:
> On Sun, Jun 9, 2013 at 7:58 PM, John Assael <[email protected]> wrote:
>
>> Dear Jonathan,
>> How are you? Your help was very valuable to me!
>> I am running a different model this time and I am using sort of a trial
>> and error to simulate the real experiment and find the values that make a
>> perfect fit to the experimental results.
>>
>> The problem is that the first 2 seconds of the simulation are the same
>> each time, as due to the mesh geometry are not affected by the different
>> parameters I use each time.
>>
>> Is there a way to save the process until that stage, saving me a lot of
>> extra time from each run?
>>
>> FiPy has no built in way to save the state of a simulation and restart
> with a single command or two. You will need to explicitly save out the
> arrays for each of the variables and then reassign the values when the
> simulation is restarted. There are lots of ways to save arrays, the
> simplest is to use "np.savetxt" and "np.loadtxt" to get started. You can
> also use "fipy.dump" and "fipy.load" to pickle the variables directly.
> PyTables is also a possibility for more sophisticated data across multiple
> time steps in one file with quick access.
>
> For the actually reading and writing, one way is to have two different
> setups. One setup reads in whatever the latest variable values are from the
> saved files while the other reinitializes. Something like this
>
> mesh = ...
>
> if restart:
>
> timestep = np.loadtxt('timestep.txt')
> var0array = np.loadtxt('var0.txt')
> var1array = np.loadtxt('var1.txt')
>
> else:
>
> timestep = 0
> var0array = np.zeros(...)
> var1array = np.zeros(...)
>
> var0 = fp.CellVariable(mesh=mesh, value=var0array)
>
> ...
>
>
> while timestep < timesteps
> np.savetxt('timestep.txt', np.array((timestep,)))
> np.savetxt('var0.txt', var0.value)
> np.savetxt('var1.txt', var1.value)
>
> ## solve equations
>
>
>>
>> Finally, the timeStep = 2e-4 works perfect but its too slow, is there a
>> better way to increase the time step, for example gradually?
>>
>
> The time step is entirely under your control you can increase it based on
> elapsed time or something more sophisticated such as the number of steps to
> convergence, but FiPy doesn't have any build in stepper schemes.
>
> --
> Daniel Wheeler
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
--
..:: ic3man.gr ::..
software | design | development
_/_/_/ email: [email protected]
_/_/_/ www : http://www.ic3man.gr
------------------------------------------------------------------
This e-mail and any attachments are confidential. You may not copy
or disseminate any information contained in them to anyone other
than the intended recipient. If you are not the intended recipient
please contact the sender by reply e-mail and destroy all copies
of the original message immediately.
------------------------------------------------------------------
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
