/*===========================================================================*/
Summary
ISAC is a static and dynamic thermal analysis software package.  Its initial
version was described in the following research paper:
  Y. Yang, Z. Gu, C. Zhu, L. Shang, and R. P. Dick, "Adaptive Chip-Package
  Thermal Analysis for Synthesis and Design," Proc. Conf. on Design,
  Automation, and Test in Europe, March 2006.
The goal of ISAC is to support rapid and accurate steady-state and dynamic
thermal analysis for use in integrated circuit synthesis and computer
architecture research.
/*===========================================================================*/
Simplest example for static thermal analysis here 
/* The setup for the example is as follows. 
The chip has one silicon layer and one heatsink layer
The size of the silicon layer is 0.8E-3 meter * 0.8E-3 meter * 0.2E-3 meter 
in [x, y, z] dimension. The thermal conductivity of the silicon layer is
157 Watt/meter-Kelvin. The Specific heat capacity of the silicon layer is
1.638E6 Joule/meter^3-Kelvin.
The size of the heatsink layer is 1.6E-3 meter * 1.6E-3 meter * 0.6E-3 meter 
in [x, y, z] dimension. The thermal conductivity of the heatsink layer is
400 Watt/meter-Kelvin. The Specific heat capacity of the heatsink layer is
3.55E6 Joule/meter^3-Kelvin.
The silicon layer has four homogeneous functional units, and each of them has 
a width/length of half those of the silicon layer.  Each of the functional 
units has power consumption.
NOTE: The 0 point of the x, y axis is set to be at the center of both of the
silicon and heatsink layer, and the 0 point for the z axis is set to be at the
top of the silicon layer. Usually, the active layer is also at the top of the 
silicon layer.
The power of the functional unit at position (-0.2E-3, -0.2E-3) is 3.0 Watt;
The power of the functional unit at position (0.2E-3, -0.2E-3) is 4.0 Watt;
The power of the functional unit at position (-0.2E-3, 0.2E-3) is 1.0 Watt;
The power of the functional unit at position (0.2E-3, 0.2E-3) is 2.0 Watt. */
/* 2x2 power grid in active layer with default power consumption of 0.0 W. */
	vector<Element> power_profile = create_grid(2, 2, 0.8E-3, 0.8E-3);
/* Assign the power consumption values to the Elements. */
	power_profile[0].value = 1.0;
	power_profile[1].value = 2.0;
	power_profile[2].value = 3.0;
	power_profile[3].value = 4.0;
/* Create the chip package. */
	vector<Layer> chip_package;
/* Silicon conductivity, specific heat capacity, sizes [x, y, z], followed by
   heatsink conductivity, specific heat capacity, sizes [x, y, z]. */
	chip_package = create_chip_package(157, 1.638E6,
		0.8E-3, 0.8E-3, 0.2E-3, 400, 3.55E6, 1.6E-3, 1.6E-3, 0.6E-3);
/* Construct solver with a silicon layers, a copper heatsink and 45 degree 
   celsius ambient temperature. */
	ISAC isac(chip_package, Boundary::cu45deg);
/* Static thermal analysis solver follows. */
	vector<Element> thermal_profile = isac.solve_static(power_profile);
/* Output the profile.
   See later examples for other output methods and conversion routines. */
  for (vector<Element>::const_iterator i = thermal_profile.begin();
		i != thermal_profile.end(); ++i)
	{
		cout << *i << endl;
	}
/*===========================================================================*/
Contents
Interface
  isac.h: Thermal analysis C++ header file
  isac_c.h: Thermal analysis C header file
  grids.h: Grid basis conversion support routines C++ header file
Implementation
  isac.cct: C++ template implementation file
  grids.cct: C++ template implementation file
  libisac.a: Library with C and C++ bindings
/*===========================================================================*/
High-level description
Chip-package configuration
The ISAC(.) constructor creates an instance of a thermal analysis object.
This constructor accepts, as input, a Boundary object specifying the
chip-package boundary conditions.  The ISAC constructor also accepts a C++
standard template library (STL) container of Layers, which describes the
construction of the chip and package.  Each material Layer has a thermal
conductivity, a specific heat capacity, a width, a height, and a depth.  The
ordering of layers in the chip and package is determined by their order in the
container.
Steady-state thermal analysis
Steady-state thermal analysis is conducted via the solve_static(.) method,
which accepts an STL container of thermal elements describing the chip-package
power consumption profile.  Each element is a rectangular parallelipiped
having a width, height, depth, a position in three-dimensional space, and a
power consumption value.  This interface allows three-dimensional power
profiles to be specified.  However, if the depth and z-axis positions are not
specified, the thermal elements default to the IC active layer, making the
common case of providing a two-dimensional active layer power profile easy.
The solve_static(.) method returns an STL vector of thermal elements, each of
which has a temperature.
Dynamic thermal analysis
Dynamic thermal analysis is conducted via the init_dynamic(.) and
step_dynamic(.) methods.  init_dynamic(.) initializes the chip-package to the
steady-state thermal profile associated with power profile argument.  This
method is similar to solve_static(.).  However, it stores the initial
temperature profile as internal state to allow subsequent dynamic analysis.
Starting from the current thermal profile state, step_dynamic(.) method
determines the new thermal profile if "duration" seconds are spent with the
power profile argument.  This temperature profile is returned from the method
and stored as internal state for further time advances.
Generic mirror interface
For the sake of convenience, we have provided a generic mirror interface
allowing users to use any STL container, including a plain array, to provide
power profiles.  This interface is otherwise identical to the primary
interface.
Grid basis conversion
Specifying power and thermal profiles using unordered containers of rectangular
parallelipipeds is general.  However, some users may want to provide power
profiles, or receive thermal profiles, in homogeneous arrays.  We have provided
a number of generic, template-based, conversion routines to convert to and from
containers of elements to two-dimensional and three-dimensional homogeneous
arrays.  We decided that it would be better to support element basis conversion
via routines than to complicate the interface with ISAC.
Recall that ISAC is spatially adaptive.  The elements in the output thermal
profiles may be heterogeneous, and may not have the same positions and sizes
as the elements in the input power profiles.  The element basis conversion
routines can also be used to find the average, minimum, and maximum
temperatures encountered within three-dimensional parallelipipeds
corresponding to, e.g., processor functional units.
/*===========================================================================*/
Demonstrating how to use the library 
/* Before conducting either the static or the dynamic thermal analysis, some
   parameters have to be set as follows (Step 1 to Step 4). */
/* Step 1 */
/* Generate the 2(or 3)-dimensional power profile array in the format of
   vector<Element>, they can be either homogeneous or heterogeneous power 
   block elements) */
   
/* Homogeneous power profile is easy to be generated by calling the routine
   create_grid(.) */
/* Heterogeneous power profile has to be specified as follows. */
/* First, generate the container for the power blocks. */ 
	vector<Element> power_profile;
/* Then, assign power values to the power block elements by doing
   the following for all the power blocks */
	power_profile.push_back(Element(double size_x, double size_y, double 
		center_x, double center_y, double value));
	or 
	power_profile.push_back(Element(double size_x, double size_y, 
		double size_z, double center_x, double center_y, double center_z, 
		double value));
/* Here, 
   size_x, size_y, size_z are the according sizes in the [x, y, z] dimension
   for each power block;
   center_x, center_y, center_z are the according center positions in the  
   [x, y, z] dimension for each power block. */
/* Step 2 */
/* Define the chip package information. Now at most two layers are supported,
   usually they are one silicon layer and one heatsink layer. 
   Note: When a heatsink layer is defined, its sizes in [x, y] dimension CANNOT 
   be smaller than those of the silicon layer. */
	vector<Layer> chip_package;
	chip_package = create_chip_package(double sil_cond, double sil_cap,
		double sil_sizex, double sil_sizey, double sil_sizez,
		double sink_cond, double sink_cap, double sink_sizex,
		double sink_sizey, double sink_sizez);
/* In the create_chip_package function, the units are:
   W-watt, m-meter, K-Kelvin temperature,
   J-Joule, Kg-Kilogram.
   The parameters are:
   sil_con-the conductivity of the silicon layer in W/m-K;
   sil_cap-the specific heat capacity of the silicon layer in J/m^3-K;
   sil_cap can be obtained by
     (heat capacity per Kg (J/Kg-K) * density (Kg/m^3)) of silicon;
   sil_sizex-the size of the silicon layer in x dimension;
   sil_sizey-the size of the silicon layer in y dimension;
   sil_sizez-the size of the silicon layer in z dimension;
   sink_con-the conductivity of the heatsink layer in W/m-K;
   sink_cap-the specific heat capacity of the heatsink layer in J/m^3-K;
   sink_cap can be obtained by
     (heat capacity per Kg (J/Kg-K) * density(Kg/m^3)) of heatsink;
   sink_sizex-the size of the heatsink layer in x dimension;
   sink_sizey-the size of the heatsink layer in y dimension;
   sink_sizez-the size of the heatsink layer in z dimension.
   If there is no heatsink layer, then only the parameters for the silicon
   layer need to be specified. The parameters of the heatsink layer will
   have default zero values under such circumstances. */
/* Step 3 (optional) */
/* Define the boundary conditions, or choosing the existing preset boundary
   conditions */
	Boundary boundary_cond(double a_t, double a_c);
   /* In the boundary_cond(.) function:
   a_t-the ambient temperature value;
   a_c-the ambient convection resistance in K/m^2-W, and another part of 
   the ambient temperature will be obtained by (a_c * heatsink area * total power).
   Therefore, the actual ambient temperature is 
   (a_t + a_c * heatsink area * total power). */
/* Step 4 */
/* Generate an object of ISAC based on the chip_package and the boundary
   conditions, to facilitate the following thermal analysis. */
	ISAC(const vector<Layer> & chip_package,
		const Boundary & boundary_cond);
/* In the above function, chip_package, and boundary_cond are generated
   in Step 2 and Step 3 respectively. */
/* Step 5 */
/* Now static, and/or dynamic thermal simulation can be conducted 
   respectively.*/
/* Step 5.a */
/* Static thermal analysis performs the following task:
   Given the power profile, to determine the stable temperature profile that
   the integrated chip will converge to when time approaches infinity. */
/* The static thermal analysis should be conducted by calling 
   the solve_static(.) function. */
	vector<Element> thermal_profile = isac.solve_static(
		const vector<Element> & power_profile);
/* Note: The power_profile has already been generated in Step 1. */
/* Or Step 5.b */
/* Dynamic thermal analysis performs the following task:
   Given the initial power profile, call the init_dynamic function to obtain
   the initial temperature profile, then call the step_dynamic function to
   determine the intermediate temperature profile at any interested duration 
   time based on the new power profile. */
/* Step 5.b.1 */
/* Init_dynamic(.) solver should be called to determine the initial temperature
   profile based on the initial power profile. */
	thermal_profile = isac.init_dynamic(
		const vector<Element> & power_profile);
/* The power_profile has already been generated in Step 1. */
/* Step 5.b.2 */
/* Define the simulation time for the dynamic thermal analysis based on
   new power profile and the existing initial temperature profile. */
	double duration = *** FILL IN VALUES***;
/* Step 5.b.3 */
/* Call the step_dynamic(.) function to do the dynamic thermal analysis, to
   determine the temperature profile of the chip after the duration time. */
	thermal_profile = isac.step_dynamic(const vector<Element> &
		power_profile, double duration);
/* Step 5.b.** */
/* Whenever a new power profile is obtained, step 5.b.3 can be repeated
   to do the dynamic thermal analysis based on the temperature
   profile obtained from the last thermal analysis step. */
/* Step 6 */
/* Output temperature values. */
/* Step 6.a */
/* After Step 5 (either Step 5.a or Step 5.b), the data structure
   thermal_profile (a vector of Element) already contains the simulated
   temperature values. Look at the 'value' property of each Element for
   the temperature values of the whole chip. */
/* Step 6.b */
/* If the temperature values of the active layer are desired, call the
   following convert_elts_to_homog_active(.) function to obtain them 
   (for homogeneous temperature output). */
/* Generate a vector of double type to store the output temperature of the
   active layer. */
	vector<vector<double> > active_temp;
/* Convert the temperature values from the thermal blocks to the above
   vector. */
	active_temp = convert_elts_to_homog_active(
		const Element * thermal_profile.begin(),
		const Element * thermal_profile.end(), int divx, int divy, 
		double sizex, double sizey);
/* In the above convert function:
   The parameter thermal_profile is obtained from Step 5;
   divx/divy means the number of blocks desired for the temperature output
   in [x,y] dimension;
   sizex/sizey means the size of the chip part desired for the temperature
   output in [x,y] dimension. */
/* Output the temperature of the desired blocks. */
	for (unsigned int i = 0; i < active_temp.size(); ++i){
		for (unsigned int j = 0; j < active_temp[i].size(); ++j){
			cout << active_temp[i][j] << " ";
		} cout << endl;
	} 
/* Step 6.c */
/* Convert the temperature of the thermal blocks into that of the elements
   in the active layer(to heterogeneous temperature output). */
/* Generate a vector of type 'Element' to store the parameters of the
   block elements desired for the temperature output, and also define
   the parameters (sizes, centers) of the desired blocks. */
	vector<Element > active_elts_temp;
	active_elts_temp.push_back(Element(double size_x, double size_y, 
		double center_x, double center_y, double value));
	or 
	active_elts_temp.push_back(Element(double size_x, double size_y, 
		double size_z, double center_x, double center_y, double center_z, 
		double value));
/* To calculate the temperature of the desired blocks based on the
   thermal_profile from the thermal analysis. */
	convert_elts_to_elts_avg(const Element * thermal_profile.begin(), 
		const Element * thermal_profile.end(),
		Element * active_elts_temp.begin(), Element * active_elts_temp.end());
/* Output the temperature of the desired blocks. */
	for (vector<Element>::const_iterator i = active_elts_temp.begin();
		i != active_elts_temp.end(); ++i)
	{
		cout << *i << endl;
	}
% Bill of goods
/* The proposed techniques make both dynamic and static thermal analysis practical
  within the inner loop of IC synthesis algorithms. */
% Summary of interface
The main interface to be taken care of: 
1. The chip_package, which contains the silicon layer and/or the heatsink
   layer.
2. The class isac, its object will be used to do the thermal analysis.
3. The Element structure. It contains the sizes information, center
   position information and the value (either power numbers or temperature
   numbers). The input power profile for the whole chip must be in the
   format of vector<Element> in order to be used in the thermal analysis
   functions. 
% Example code for each mode of use
/*===========================================================================*/
Example code
/*###########################################################################*/
/* Example code for static thermal analysis */ 
/* The chip has two layers, one silicon layer and one heatsink layer
The size of the silicon layer is 0.8E-3 meter * 0.8E-3 meter * 0.2E-3 meter in
[x, y, z] dimension. The thermal conductivity of the silicon layer is
157 Watt/meter-Kelvin. The Specific heat capacity of the silicon layer is
1.638E6 Joule/meter^3-Kelvin.
The size of the heatsink layer is 1.6E-3 meter * 1.6E-3 meter * 0.6E-3 meter in
[x, y, z] dimension. The thermal conductivity of the heatsink layer is
400 Watt/meter-Kelvin. The Specific heat capacity of the heatsink layer is
3.55E6 Joule/meter^3-Kelvin.
NOTE: The 0 point of the x/y axis is set to be at the center of the
silicon chip.
Heterogeneous power profile as follows.
The power of the functional unit (size 0.4E-3*0.4E-3) at position 
(-0.2E-3, -0.2E-3) is 1.0 Watt;
The power of the functional unit (size 0.4E-3*0.4E-3) at position 
(0.2E-3, -0.2E-3) is 2.0 Watt;
The power of the functional unit (size 0.4E-3*0.4E-3) at position 
(-0.2E-3, 0.2E-3) is 3.0 Watt;
The power of the functional unit (size 0.2E-3*0.4E-3) at position 
(0.1E-3, 0.2E-3) is 2.0 Watt;
The power of the functional unit (size 0.2E-3*0.4E-3) at position 
(0.3E-3, 0.2E-3) is 2.0 Watt. */
/* Step 1 */
/* Generate the power profile container. */
	vector<Element> power_profile;
/* Assign sizes, centers, and power consumption numbers for each element. */
	power_profile.push_back(Element(0.4E-3, 0.4E-3, -0.2E-3, -0.2E-3, 1.0));
	power_profile.push_back(Element(0.4E-3, 0.4E-3, 0.2E-3, -0.2E-3, 2.0));
	power_profile.push_back(Element(0.4E-3, 0.4E-3, -0.2E-3, 0.2E-3, 3.0));
	power_profile.push_back(Element(0.2E-3, 0.4E-3, 0.1E-3, 0.2E-3, 2.0));
	power_profile.push_back(Element(0.2E-3, 0.4E-3, 0.3E-3, 0.2E-3, 2.0));
/* Step 2 */
/* Create the chip package. */ 
	vector<Layer> chip_package;
	chip_package = create_chip_package(157, 1.638E6,
		0.8E-3, 0.8E-3, 0.2E-3, 400, 3.55E6, 1.6E-3, 1.6E-3, 0.6E-3);
/* Step 3 is omitted by using preset boundary conditions. */
/* Step 4 */
/* Construct solver with a copper heatsink and 45 degree celsius ambient
   temperature. */
	ISAC isac(chip_package, Boundary::cu45deg);
/* Step 5 */
/* Conduct the static thermal analysis. */
	vector<Element> thermal_profile = isac.solve_static(power_profile);
/* Step 6.a */
/* Output the temperature of the whole chip blocks. */
	for (vector<Element>::const_iterator i = thermal_profile.begin();
		i != thermal_profile.end(); ++i)
	{
		cout << *i << endl;
	}
/* Step 6.b */
/* If the temperature value of the active layer is desired, call the
   following convert_elts_to_homog_active(.) function to obtain them 
   (for homogeneous temperature output). */
/* We are going to output the temperature values for the active layer with 
   16 homogeneous temperature grid */
/* Generate a vector of double type to store the output temperature of the
    active layer. */
	vector<vector<double> > active_temp;
/* Convert the temperature values from the thermal blocks to the above 
   vector. */
	active_temp = convert_elts_to_homog_active(thermal_profile.begin(),
		thermal_profile.end(), 4, 4, 0.8E-3, 0.8E-3);
/* Output the temperature values in 4*4 blocks. */
		for (unsigned int i = 0; i < active_temp.size(); ++i){
			for (unsigned int j = 0; j < active_temp[i].size(); ++j){
				cout << active_temp[i][j] << " ";
			}	cout << endl;
		}
/* Step 6.c */
/* To output the temperature of each power block in the active layer. */
/* First generature the tempearture container of the power blocks. */
	vector<Element> active_elts_temp;
	for (int i = 0; i < power_profile.size(); ++i){
		active_elts_temp.push_back(power_profile[i]);
		active_elts_temp[i].value = 0.0;
	}
	convert_elts_to_elts_avg(thermal_profile.begin(),
		thermal_profile.end(), active_elts_temp.begin(), active_elts_temp.end());
/* Output the temperature of the active blocks. */
	for (vector<Element>::const_iterator i = active_elts_temp.begin();
		i != active_elts_temp.end(); ++i)
	{
		cout << *i << endl;
	}
/*###########################################################################*/
/* Example code for dynamic thermal analysis */
/* The chip only has two layers: one silicon layer and one heatsink layer
The size of the silicon layer is 0.8E-3 meter * 0.8E-3 meter * 0.2E-3 meter in
x/y/z dimension. The thermal conductivity of the silicon layer is
157 Watt/meter-Kelvin. The Specific heat capacity of the silicon layer is
1.638E6 Joule/meter^3-Kelvin.
The size of the heatsink layer is 1.6E-3 meter * 1.6E-3 meter * 0.6E-3 meter in
[x y z] dimension. The thermal conductivity of the heatsink layer is
400 Watt/meter-Kelvin. The Specific heat capacity of the heatsink layer is
3.55E6 Joule/meter^3-Kelvin.
The silicon layer has four homogeneous functional units, and each of them has a
width/length of half those of the silicon layer.  Each of them has a
power number.
NOTE: The 0 point of the x/y axis is set to be at the center of the
silicon chip.
Initial heterogeneous power consumption:
The power of the functional unit (size 0.4E-3*0.4E-3) at position 
(-0.2E-3, -0.2E-3) is 1.0 Watt;
The power of the functional unit (size 0.4E-3*0.4E-3) at position 
(0.2E-3, -0.2E-3) is 2.0 Watt;
The power of the functional unit (size 0.4E-3*0.4E-3) at position 
(-0.2E-3, 0.2E-3) is 3.0 Watt;
The power of the functional unit (size 0.2E-3*0.4E-3) at position 
(0.1E-3, 0.2E-3) is 2.0 Watt;
The power of the functional unit (size 0.2E-3*0.4E-3) at position 
(0.3E-3, 0.2E-3) is 2.0 Watt.
The power profile has been used to call the init_dynamic(.) function to
obtain the initial temperature profile. After that, the power of each FU 
increases by 0.1 Watt, and the step_static(.) function is called to 
calculate the transient temperature profile at the time slot specified 
by parameter 'duration'. */
/* Step 1 */
/* Generate the power profile container. */
	vector<Element> power_profile;
/* Assign the sizes, centers and the power numbers for each Element */
	power_profile.push_back(Element(0.4E-3, 0.4E-3, -0.2E-3, -0.2E-3, 1.0));
	power_profile.push_back(Element(0.4E-3, 0.4E-3, 0.2E-3, -0.2E-3, 2.0));
	power_profile.push_back(Element(0.4E-3, 0.4E-3, -0.2E-3, 0.2E-3, 3.0));
	power_profile.push_back(Element(0.2E-3, 0.4E-3, 0.1E-3, 0.2E-3, 2.0));
	power_profile.push_back(Element(0.2E-3, 0.4E-3, 0.3E-3, 0.2E-3, 2.0));
/* Step 2 */
/* Create the chip package. */
	vector<Layer> chip_package;
	chip_package = create_chip_package(157, 1.638E6,
		0.8E-3, 0.8E-3, 0.2E-3, 400, 3.55E6, 1.6E-3, 1.6E-3, 0.6E-3);
/* Step 3 is omitted by using default boundary conditions */
/* Step 4 */
/* Construct solver with a copper heatsink and 45 degree celsius ambient
   temperature. */
	ISAC isac(chip_package, Boundary::cu45deg);
/* Step 5 */
   /* Do the initial thermal analysis to obtain the initial temperature 
   profile. */
	vector<Element> thermal_profile = isac.init_dynamic(power_profile);
    
   /* Generate power variation here. */
	for (unsigned int i = 0; i < power_profile.size(); ++i){
		power_profile[i].value += 0.1;
	}
   /* Define the duration time desired for the dynamic thermal simulation. */
	double duration = 1E-8;
	thermal_profile = isac.step_dynamic(power_profile, duration);
   /* Whenever a new power profile is obtained, step_dynamic(.) can be called
   to do the dynamic thermal analysis based on the temperature
   profile obtained from the last thermal analysis step. */
 
/* For Step 6.a, 6.b and 6.c, please refer to those in the above static 
   analysis example. */