Remeshing and code speedup
Hi, Thanks for the help with getting FiPy running under Linux! I am trying to re-create a 1D nonlinear diffusion problem for which we have C++ code that uses the implicit Thomas algorithm based on J. Weickert, B. Romerny, M. Viergever, "Efficient and Reliable Schemes for Nonlinear Diffusion Filtering”, IEEE transactions on Image Processing, vol.7, N03, page 398, March 1998 I have been able to get results in FiPy that match this code very closely which was a great start. Our C++ code uses a fixed number of spatial points and a fixed time step, but re-meshes space to most efficiently use the size of the array; it increases the spatial step size by 2 whenever the concentration at a particular point reaches a set threshold. I tried implementing this in FiPy as well, but haven’t had much luck so far. I saw an old mailing-list entry from 2011 where a user was told that FiPy wasn’t meant to do remeshing. Is that still the case? I’d imagine one would somehow need to update the Grid1D object with the new ‘dx’, but since the CellVariable that holds the solution was initialized with that mesh object, I am not sure that such a change would propagate in a sensible fashion. I think I know how to map the value of the CellVariable to account for the change in ‘dx’ by array_size = 2000 phi.value = numpy.concatenate((phi.value[1:array_size/2:2], numpy.zeros(1500))) for the case when the initial variable holds 2000 spatial points. Maybe there’s a more elegant way, but I think this works in principle. Another question would be execution speed. Right now, even when not plotting the intermediate solutions, it takes many seconds on a very powerful computer to run a simple diffusion problem. I am probably doing something really wrong. I wasn’t expecting the code to perform as well as the C++ code, but I had hoped to come within an order of magnitude. Are there ways to optimize the performance? Maybe select a particularly clever solver? If someone could point me into the right direction that’d be great. In the end, I would like to expand the code into 2D, but given the poor 1D performance, I don’t think that this would be feasible at this point. Thanks, Carsten _Dipl.-Phys. Carsten Langrock, Ph.D. Senior Research Scientist Edward L. Ginzton Laboratory, Rm. 202 Stanford University 348 Via Pueblo Mall 94305 Stanford, CA Tel. (650) 723-0464 Fax (650) 723-2666 Ginzton Lab Shipping Address: James and Anna Marie Spilker Engineering and Applied Sciences Building 04-040 348 Via Pueblo Mall 94305 Stanford, CA _ ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: thermal contact between two regions
On Mon, Jul 23, 2018 at 1:06 PM, Drew Davidson wrote: > > 'dx' is actually the thickness of the thermal contact region. It is not the > volume of the cell. If it were a volume, the ratio dx/mesh._cellDistances > would not be dimensionless. I was seeking a thermal contact region of zero > thickness (dx=0), giving a true thermal contact resistance. Sorry, my mistake on that one. > I let this thread drop because the Salazar paper I referenced gave me > expressions for effective thermal conductivity and effective thermal > diffusivity, and this allowed me to arrive at a reasonable-looking FiPy > result for transient thermal behavior. My result still deviates a bit from > the exact solution at steady state, and this deviation persists despite mesh > refinement, but I am busy with other things and have called it good enough > for now. Thanks for following up and sharing your code. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: thermal contact between two regions
Hello Dr. Wheeler,
'dx' is actually the thickness of the thermal contact region. It is not
the volume of the cell. If it were a volume, the ratio
dx/mesh._cellDistances would not be dimensionless. I was seeking a thermal
contact region of zero thickness (dx=0), giving a true thermal contact
resistance.
I let this thread drop because the Salazar paper I referenced gave me
expressions for effective thermal conductivity and effective thermal
diffusivity, and this allowed me to arrive at a reasonable-looking FiPy
result for transient thermal behavior. My result still deviates a bit from
the exact solution at steady state, and this deviation persists despite
mesh refinement, but I am busy with other things and have called it good
enough for now.
OK, I won't leave you hanging. Here is the code I was working on. Maybe
you or someone else can spot a defect that would explain why the code does
not exactly replicate the steady state solution. Another thing to do is:
compare the transient result with the transient result obtained another way
(ANSYS etc). I used to have an analytical solution for the transient
laying around, but threw it out in a pile of papers (oops).
#I am trying to look at the script by Wheeler at
#https://www.mail-archive.com/fipy@nist.gov/msg02640.html
#and modify it in view of Guyer as in
#https://www.mail-archive.com/fipy@nist.gov/msg02641.html
#now instead of a steady state problem, do a transient problem
import fipy as fp
import numpy as np
from IPython import embed
import sys
import os
import pygmsh
import meshio
import csv
import matplotlib.pyplot as plt
from shutil import copyfile
#head
def
get_uniform_mesh_with_pygmsh(L_Sample,L_Substrate,cellSize,write_to_msh_file=False):
'''
figure it out in ipython
ref: https://github.com/nschloe/pygmsh/blob/master/test/test_pacman.py
'''
geom = pygmsh.built_in.Geometry()
#geom.add.point? see what the input arguments are in ipython using ?
p1=geom.add_point([0,0,0],cellSize)
p2=geom.add_point([L_Substrate,0,0],cellSize)
p3=geom.add_point([L_Substrate+L_Sample,0,0],cellSize)
l1=geom.add_line(p1,p2)
l2=geom.add_line(p2,p3)
#not much point in doing physical groups since have no way to transmit
this info to 1D fipy mesh
# geom.add_physical_point(p1,label='SubstrateEnd')
# geom.add_physical_point(p2,label='TCR')
# geom.add_physical_point(p3,label='SampleEnd')
# geom.add_physical_line(l1,label='Substrate')
# geom.add_physical_line(l2,label='Sample')
points, cells, point_data, cell_data, field_data =
pygmsh.generate_mesh(geom)
if write_to_msh_file:
mesh=meshio.Mesh(points=points,cells=cells,point_data=point_data,cell_data=cell_data,field_data=field_data)
meshio.write('TCR01.msh',mesh)
mesh=get_fipy_1D_mesh_from_points(points)
return mesh
#head
#head
if __name__=='__main__':
#substrate is copper and sample is pure tin
rho_Sample=6980. #use wolfram alpha to get the density of copper;
kg/m^3
cp_Sample=227. #use wolfram alpha to get the specific heat of copper;
J/kg/K
k_Sample=59.6 #W/m/K
D_Sample=k_Sample/rho_Sample/cp_Sample
L_Sample=2e-6 #m SETTING
rho_Substrate=8960. #use wolfram alpha to get the density of copper;
kg/m^3
cp_Substrate=384.4 #use wolfram alpha to get the specific heat of
copper; J/kg/K
k_Substrate=381.8 #W/m/K
D_Substrate=k_Substrate/rho_Substrate/cp_Substrate
L_Substrate=10*L_Sample #m SETTING
#model an interface thermal resistance between copper and solder, with
the solder being soldered to the copper
R_TIM_unit_area=2e-6 #K/(W/m2) SETTING
TIM_target_dT=1. #K SETTING
HeatFlux=TIM_target_dT/R_TIM_unit_area #W/m^2 heat flux is determined
by TIM_target_dT
#using a uniform mesh:
cellSize=.1*L_Sample #SETTING
mesh=get_uniform_mesh_with_pygmsh(L_Sample,L_Substrate,cellSize,write_to_msh_file=False)
#make sure you like the mesh before doing a run
#sys.exit()
diffCoeff = fp.FaceVariable(mesh=mesh, value=D_Substrate)
diffCoeff.setValue(D_Sample,where=(mesh.faceCenters[0]>=L_Substrate))
xCenters,=mesh.cellCenters
#thermal conductivity of a cell if it's R" were R_TIM_unit_area:
Kcontact=1./R_TIM_unit_area*mesh._cellDistances #type is
numpy.ndarray; [W/m/K]
Kcell=fp.CellVariable(mesh=mesh,value=k_Substrate)
Kcell.setValue(k_Sample,where=(xCenters>=L_Substrate))
Kavg = Kcell.harmonicFaceValue #a FaceVariable
#see 20180717EffectiveProperties.xoj
Keff=Kcontact * Kavg / (Kavg + Kcontact) #W/m/K
rho_times_cp_eff=.5*(rho_Substrate*cp_Substrate+rho_Sample+cp_Sample)
#TODO ask on fipy mailing list as to how to express dAP1/(dAP1+dAP2) in
fipy syntax; the factor of .5 here is only approximate
diffCoeff.setValue(Keff/rho_times_cp_eff, where=(mesh.faceCenters[0] ==
L_Substrate))
v = fp.CellVariable(mesh=mesh)
v.constrain(0, where=mesh.facesLeft)
v.faceGrad.constrain(HeatFlux /
Re: thermal contact between two regions
On Tue, Jul 17, 2018 at 4:43 PM, Drew Davidson wrote: > > I still need FiPy code for: > > > dAP1/(dAP1+dAP2) > > > and > > > dAP2/(dAP1+dAP2) > > > where dAP1 and dAP2 were distances from cell center to cell face for cells > on either side of the interface. If these FiPy expressions are unavailable, > I would think assuming .5 is going to be OK... I think it's outlined in the link that you gave. Kcontact = hc * mesh._cellDistances Kavg = Kcell.harmonicFaceValue K.setValue(Kcontact * Kavg / (Kavg + Kcontact * (1 - dx / mesh._cellDistances)), where=mesh.physicalFaces['thermal contact']) m._cellDistances is the distance between each cell and `dx` is the volume of the cell for a 1D problem. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
