dropbox is completely blocked at NIST. I've sent you a link to enable file upload.
On Sep 9, 2014, at 1:19 PM, John Assael <[email protected]> wrote: > 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 ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
