I have uploaded all the files at:
https://www.dropbox.com/sh/1wtpic28i8p8o45/AAARJRJn05X4VbvESJrGqL_1a?dl=0
fipy_tmp.py - is the main script that produces the difference of
temperature at the radius of "mat1"
mesh_com.msh - the mesh from comsol with average element quality 94%
expected_results.txt - and the expected results that I should be getting
It is worth mentioning that in my case the time steps differ a little so
check the ones every 0.05sec that match.
Thank you very much for your help it is deeply appreciated!
a copy of the script can be found below:
-------------------------------
# START Global variables
time_steps = np.array(
[3.6013E-06, 3.6013E-06, 7.2025E-06, 1.4405E-05, 2.8810E-05,
2.8810E-05, 2.8810E-05, 5.7620E-05, 5.1858E-05,
5.1858E-05, 5.1858E-05, 1.0372E-04, 1.0372E-04, 2.0743E-04,
2.0743E-04, 2.0743E-04, 2.0743E-04, 2.0743E-04,
2.0743E-04, 4.1487E-04, 4.1487E-04, 4.1487E-04, 4.1487E-04,
4.1487E-04, 8.2973E-04, 8.2973E-04, 8.2973E-04,
8.2973E-04, 8.2973E-04, 1.6595E-03, 1.6595E-03, 1.6595E-03,
1.6595E-03, 1.6595E-03, 1.6595E-03, 3.3189E-03,
3.3189E-03, 3.3189E-03, 3.3189E-03, 3.3189E-03, 6.6378E-03,
6.6378E-03, 2.1695E-03, 4.3390E-03, 4.3390E-03,
4.3390E-03, 8.6781E-03, 8.6781E-03, 8.6781E-03, 1.0949E-02,
1.0949E-02, 2.1897E-02, 1.7154E-02, 2.5076E-02,
2.4924E-02, 2.4924E-02, 2.4924E-02, 1.5225E-04, 3.0450E-04,
6.0900E-04, 1.2180E-03, 2.4360E-03, 4.8720E-03,
9.7440E-03, 1.9488E-02, 1.1328E-02, 2.2657E-02, 2.7343E-02, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05,
0.05, 0.05])
time_steps_cum = np.cumsum(time_steps)
# END Global variables
def calcDT(mesh, properties, debug=False):
global s_dT
global s_values_time
global s_values_timeLog
global s_values_dt
t0 = time.time()
values_time = []
values_dt = []
x, y = mesh.cellCenters
x_face, y_face = mesh.faceCenters
if debug:
print 'Start Init Variables', time.time() - t0
dT = CellVariable(name=r'$\Delta T$', mesh=mesh, value=0.)
QCell = CellVariable(name=r'$Q$', mesh=mesh, value=0.)
tcCell = FaceVariable(name=r'$Tc$', mesh=mesh, value=0.)
dencpCell = CellVariable(name=r'$DenCp$', mesh=mesh, value=0.)
# Set Mat1 Properties
meshMat1 = (x ** 2 + y ** 2 <= properties['mat1']['r'] ** 2)
meshMat1Face = (x_face ** 2 + y_face ** 2 <= properties['mat1']['r'] **
2)
QCell.setValue(properties['mat1']['Q'], where=meshMat1)
tcCell.setValue(properties['mat1']['tc'], where=meshMat1Face)
dencpCell.setValue(properties['mat1']['den'] *
properties['mat1']['cp'], where=meshMat1)
# Set Pt Properties
meshSolidPt = (y <= (properties['pt']['thick'] / 2)) & (
y >= -(properties['pt']['thick'] / 2)) & ~meshMat1
meshSolidPtFace = (y_face <= (properties['pt']['thick'] / 2)) & (
y_face >= -(properties['pt']['thick'] / 2)) & ~meshMat1Face
tcCell.setValue(properties['pt']['tc'], where=meshSolidPtFace)
dencpCell.setValue(properties['pt']['den'] * properties['pt']['cp'],
where=meshSolidPt)
# Set Mat2 Properties
meshSolidMat2 = (y > (properties['pt']['thick'] / 2))
meshSolidMat2Face = (y_face > (properties['pt']['thick'] / 2))
tcCell.setValue(properties['mat2']['tc'], where=meshSolidMat2Face)
dencpCell.setValue(properties['mat2']['den'] *
properties['mat2']['cp'], where=meshSolidMat2)
# Set Int Properties
meshSolidInt = (y < -(properties['pt']['thick'] / 2)) & (
y >= -(properties['pt']['thick'] / 2 + properties['int']['thick']))
meshSolidIntFace = (y_face < -(properties['pt']['thick'] / 2)) & (
y_face >= -(properties['pt']['thick'] / 2 +
properties['int']['thick']))
tcCell.setValue(properties['int']['tc'], where=meshSolidIntFace)
dencpCell.setValue(properties['int']['den'] * properties['int']['cp'],
where=meshSolidInt)
# Set Measure Properties
meshSolidMeasure = (y < -(properties['pt']['thick'] / 2 +
properties['int']['thick']))
meshSolidMeasureFace = (y_face < -(properties['pt']['thick'] / 2 +
properties['int']['thick']))
tcCell.setValue(properties['measure']['tc'], where=meshSolidMeasureFace)
dencpCell.setValue(properties['measure']['den'] *
properties['measure']['cp'], where=meshSolidMeasure)
# Heat Equation
heat_eq = (TransientTerm(coeff=dencpCell) ==
DiffusionTerm(coeff=tcCell) + QCell)
elapsed_time = 0.
j = 1
mySolver = LinearLUSolver()
while round(elapsed_time, 4) < 10:
time_step = time_steps[j - 1]
heat_eq.solve(var=dT, solver=mySolver, dt=time_step) #
val_dt = dT(((properties['mat1']['r'],), (0,)), order=1)[0]
elapsed_time += time_step
values_time.append(elapsed_time)
values_dt.append(val_dt)
if debug:
print '%s.\t%.10f\t%.10f\t%s' % (j, round(elapsed_time, 10),
time_step, val_dt)
j += 1
return {'time': values_time, 'dT': values_dt}
On Tue, Sep 9, 2014 at 7:32 PM, Daniel Wheeler <[email protected]>
wrote:
> On Tue, Sep 9, 2014 at 12:02 PM, John Assael <[email protected]> wrote:
>
>> Amazing! Changing T to a FaceVariable it really improved my results!
>> However, they still have a considerable gap from the optimal, I use scipy
>> LinearLUSolver, is there anything else I could do?
>>
>
> Maybe if you post your script (as simplified as possible) with the
> expected result we could look at it and see if there is anything obvious.
>
> --
> 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 ]
