Numerical methods to calculate solidification/melting fronts

[FleshJeb]

28 minutes ago, Cunjo Carl said:

Presently walls not maintained at boundary conditions, just initially set to them. Is that intended?

Was I correct about the top and bottom walls? Per her diagram at the top, they're insulated, not actively cooled or heated. Do they basically have no effect on the environment (just calc the neighbors that are NOT out of bounds), or is there some weird reflection effect that has to be taken into account?

Doh! I can't believe I missed calculating the neighbors.

[Cunjo Carl]

50 minutes ago, FleshJeb said:

Was I correct about the top and bottom walls? Per her diagram at the top, they're insulated, not actively cooled or heated. Do they basically have no effect on the environment (just calc the neighbors that are NOT out of bounds), or is there some weird reflection effect that has to be taken into account?

Good question. Fortunately no reflection for normal insulators, just no heat flow. So as you say they have no effect. Just for academic interest, atomically smooth surfaces can have heat (phonon) reflections.

[Nivee~]

5 hours ago, FleshJeb said:

Per her diagram at the top, they're insulated, not actively cooled or heated.

I am SO sorry for the misunderstanding! Actually, this is my calculation domain for the simple 2D heating problem:

Does this code make sense now?

#include<stdio.h>

int main()
{
	double T_old[100][100];
	double T_new[100][100];
    double Tin;
	int x = 0;
	int y = 0;
	printf("\nEnter the Initial Temperature: ");
	scanf("%d",Tin);
	for (x = 1; x <= 98; x++){
		for (y = 1; y <= 98; y++){
			T_old[x][y] = Tin;
			T_new[x][y] = Tin;
		}
	}

	//Boundary Conditions
    x = 0;
    for (y = 0; y <= 99; y++){
        T_old[x][y] = 85.00;
    }//LeftWall

    x = 99;
    for (y = 0; y <= 99; y++){
        T_old[x][y] = 30.00;
    }//RightWall

    y = 0;
    for (x = 0; x <= 99; x++){
        T_old[x][y] = 30.00;
    }//BottomWall
    y = 99;
    for (x = 0; x <= 99; x++){
        T_old[x][y] = 30.00;
    }//TopWall

	double alpha = 0.000137;
	double cellSize = 0.01;
	double timeStep = 0.0001;
	double coefficient = alpha*timeStep/(cellSize*cellSize);

	double time = 0.00;

	for (time = 0.00; time <= 300.00; time += timeStep){
		for (x = 1; x <= 99; x++ ){
			for (y = 1; y <= 99; y++){
				T_new[x][y] += coefficient*(T_old[x][y-1]-T_old[x][y]);
			}
		}
		for (x = 0; x <= 99; x++){
			for (y = 0; y <= 99; y++){
				T_old[x][y] = T_new[x][y];
			}
		}
	}
	for (x = 0; x <= 99; x++){
		for (y = 0; y <= 99; y++){
			printf("%d\n",T_new[x][y]);
		}
	}

	return 0;
	getch();
}

 

Edited June 5, 2019 by Nivee~

[satnet]

That code makes more sense, but you're still only calculating heat transfer for one neighbor. You have a constant heat source at x=0, which you have effectively modeled by starting at x=1 in your loop (once you have an x-1 neighbor calculation the x=1 cells will draw from the x=0 cells, but the x=0 cells will never change their value making them a constant source).  You are effectively giving y=0 a constant temperature by starting at y=1, but y=0 is still participating in heat flow so it is not an insulator. I would assume you would calculate from y=0 to y=99 and use boundary checks to just not calculate heat flow in the out of bounds direction to simulate an insulator rather than have a row of constant temperature cells.

When printing a double you should use %f or %e. For this I would use %e (scientific notation format) probably with a precision specifier (i.e. %.3e, which will give you a number with 3 numbers after the decimal like 1.234e+1). You might find %f more readable in early stages, though since it generally won't use an exponent. The %d specifier is for integers.

[FleshJeb]

@Nivee~ Basically what satnet said (and Carl on the prior page). Although if the edges are all supposed to be constant temperature, I think you can avoid testing for mathematical boundary conditions and just run your problem domain from 1 to 98 instead of 0 to 99. I think you mentioned this somewhere above.

The other thing that confuses me: Does the single alpha variable account for all methods and rates of heat transfer in the system? One of Carls early posts said:

Quote

calculate heat flow from conduction and add/subtract the enthalpy from the cells in "new"
calculate heat flow from convection and add/subtract the enthalpy from the cells in "new"
calculate heat flow from viscous friction and add/subtract the enthalpy from the cells in "new"
etc. etc, going through all of our governing equations also including for mass flow.

Programming-wise, I'd add some global variables to make changing the parameters easier.

int X_BOUND = 99; // Maximum X index. Or you could scanf.
int Y_BOUND = 99; // Maximum Y index. Or you could scanf.

// Comment your code, jerk. 

int main()
{
	double T_old[X_BOUND+1][Y_BOUND+1];
	double T_new[X_BOUND+1][Y_BOUND+1];
	double Tin;
	int x = 0;
	int y = 0;
	printf("\nEnter the Initial Temperature: ");
	scanf("%d",Tin);
	for (x = 1; x <= X_BOUND-1; x++){
		for (y = 1; y <= Y_BOUND-1; y++){

Be a good scientist:

	y = Y_BOUND;
	for (x = 0; x <= X_BOUND; x++){
		T_old[x][y] = 30.00;  // deg C
	}//TopWall

	double alpha = 0.000137;  // How was this calculated, what are the units?
	double cellSize = 0.01;   // m^2
	double timeStep = 0.0001; // seconds

Changing your output slightly will make it easier to cut and paste into a spreadsheet, or write to a CSV file.

for (x = 0; x <= X_BOUND; x++){
	for (y = 0; y <= Y_BOUND; y++){
		printf("%d,",T_new[x][y]);
	}
	printf("\n");
}

 

[Cunjo Carl]

7 hours ago, FleshJeb said:

The other thing that confuses me: Does the single alpha variable account for all methods and rates of heat transfer in the system? One of Carls early posts said:

Unfortunately everything requires its own way to be handled. Conduction is the only easy one! Rather than just having one thing we're tracking for each cell (one 'state variable'), we'll soon have several.

• Enthalpy  (rather than temp, enthalpy is a bit like the total chemical+thermal+pressure energy in the cell and is easier to use when we have different phases around)
• Pressure
• Composition (mass or mol fraction of each component)
• Velocity x,y (or momentum x,y)
• Crystal type and facing (for mush or solid cells)

 All of the above need their own equations to update them each timestep like temperature does now. To help with calculations, each timestep we'll use the above to calculate the following for each frame and also store that data:

• Temperature
• Mass fractions of liquid vs solid phase (for mush cells)
• Separately, Compositions of liquid and solid phases (for mush cells)
• Density

Finally is calculating the net force (x,y) on each cell given the pressures, density (+gravity), and shear.

The calculations for exchange between cells will always be simple linear equations like they are now, but calculating those things in the second paragraph within a cell will require simultaneous equations which are a huge slow down (I use heuristics and lookup tables instead). As a note, if we have pure Gallium, life will be much simpler! And with the assumption of one crystal type with the same conduction / solidification in all directions life becomes simpler still. @Nivee~'s paper goes through their approach on which equations to use for that case.

[Selective Genius]

@Cunjo Carl, I might be wrong, but there's something iffy about these lines in @Nivee~'s code.

for (time = 0.00; time <= 300.00; time += timeStep){
		for (x = 1; x <= 99; x++ ){
			for (y = 1; y <= 99; y++){
				T_new[x][y] += coefficient*(T_old[x][y-1]-T_old[x][y]);
			}
		}

The program calculates the temperature in a cell by using values from the neighboring cell. BUT! Only the cell at the east of the current cell. The cells on the north and south of the cell are not taken into account? Doesn't this approach reduces the 2D problem into several 1-D problems, with no connection between adjacent rows?

Something like this(Hijacking Nivee's image):

There's no connection between 1,0 and 1,1 as per as those lines of code.

Moreover,  the usual way is to discretize the domain and the equation, and then apply the TDMA method to solve the linear system of equations thus obtained, right? Or am I missing something here??

[Cunjo Carl]

16 hours ago, Selective Genius said:

@Cunjo Carl, I might be wrong, but there's something iffy about these lines in @Nivee~'s code.

...

Moreover,  the usual way is to discretize the domain and the equation, and then apply the TDMA method to solve the linear system of equations thus obtained, right? Or am I missing something here??

You're correct that the other directions are missing. It's been mentioned, so they're probably plugging away at it. Thanks for noticing and nice picture hack!

This overall shape of the FEA was the dummies version I learned back in chemical engineering undergrad, sans some pieces I didn't talk about about for calculating accumulated errors. Besides a couple odd balls, it's actually all I've ever used for my own FEAs, and though easy to apply and debug I could imagine it's not the fastest. I didn't consider matrix solution approaches, would they be appropriate here? If so I'd love to learn!

[Nivee~]

The program is being very weird.... I ran the following code:

Spoiler

#include<stdio.h>

int main()
{
	double tin;
//	double time = 10;
	double t_new[10][10];
	double t_old[10][10],t_inter[10][10];
	for (int x = 1;x < 10; x++)
	{
		for (int y = 1;y < 10; y++)
		{
                t_old[x][y]=30; //initial temperature of the interior nodes
                t_new[x][y]=30; //initial temperature of the interior nodes
		}
	}
    int x = 0;
    int y;
    for (y = 0; y <= 9; y++){
        t_old[x][y] = 30.00;
    }//Left Wall

    x = 9;
    for (y = 0; y <= 9; y++){
        t_old[x][y] = 30.00;
    }//Right Wall

    y = 0;
    for (x = 0; x <= 9; x++){
        t_old[x][y] = 85.00;
    }//Top Wall


    y = 9;
    for (x = 1; x <= 9; x++){
        t_old[x][y] = 30.00;
    }//Bottom Wall

    printf("\n---------Initial Conditions----------\n");


    for (x = 0; x< 10; x++)
	{
		for (y = 0; y < 10; y++)
		{
			printf("%d %d %f\t",x,y,t_old[x][y]);
		//	printf("    %f",t_new[x][y]);
		}
		printf("\n");
	}

	printf("\n\n\n\n\n\n");
	printf("\n---------Final Conditions----------\n");
  	float alpha = 0.000137; //Diffusivity
	float cellSize = 0.1;	//Cell size
	float timeStep = 0.01;	//
	double coefficient = alpha*timeStep/(cellSize*cellSize);

	double t = 0.00; // counter for time

	for (t = 0.00; t <= 300.00; t+= timeStep){
		for (x = 1; x <= 9; x++ ){
			for (y = 1; y <= 9; y++){
				t_new[x][y] += coefficient*(t_old[x][y-1]-t_old[x][y]);
			}
		}
		for (x = 0; x <= 9; x++){
			for (y = 0; y <= 9; y++){
				t_old[x][y] = t_new[x][y];
			}
		}
	}

	for (x = 0; x< 10; x++)
	{
		for (y = 0; y < 10; y++)
		{
			printf("%f\t",t_old[x][y]);
		}
		printf("\n");
	}


}

 

And this is what I got:

Spoiler

0.000000        0.000000        0.000000        0.000000        0.000000        0.000000        0.000000        0.000000        0.000000        0.000000
0.000000        0.492218        2.515390        6.673183        12.369412       18.222146       23.032816       26.327826       28.262227       29.255869
0.000000        0.492218        2.515390        6.673183        12.369412       18.222146       23.032816       26.327826       28.262227       29.255869
0.000000        0.492218        2.515390        6.673183        12.369412       18.222146       23.032816       26.327826       28.262227       29.255869
20238335338031785000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000       19906362662096847000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000       18541768131066443000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
15737235863061561000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000       11894756711612313000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000       7946463619192883400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
4700955072585959400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000        2477847871180646500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000        1172648712671998400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
502167307032158550000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
0.000000        0.492218        2.515390        6.673183        12.369412       18.222146       23.032816       26.327826       28.262227       29.255869
0.000000        0.492218        2.515390        6.673183        12.369412       18.222146       23.032816       26.327826       28.262227       29.255869
0.000000        0.492218        2.515390        6.673183        12.369412       18.222146       23.032816       26.327826       28.262227       29.255869
0.000000        0.492218        2.515390        6.673183        12.369412       18.222146       23.032816       26.327826       28.262227       29.255869
9581761052956567300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000 9424589882314845000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000 8778527915650223900000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000 7450733034871149200000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000 5631526244134483300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
3762222255109793200000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000 2225648872491860300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000 1173127425288470900000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000 555185966444850560000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000  237749155956853250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000

 

It looks like the top and bottom rows are disconnected, as @Selective Genius said...

[Cunjo Carl]

There's been enough back and forth that I figured I should make a little example FEA to demonstrate the idea. And just to reiterate from before, this way of doing it is quick and simple but computationally silly! Still, excepting for some error accumulation (which can be estimated) and a few oddities (like not all directions being the same in some specific cases), it works well enough to look good which was kinda the goal!

www.paikogame.com/simplefea/SimplestFEA.html

I did it in Javascript for portability- the actual FEA part should be close enough to C++ that the ideas can come across. In the demo here, you can click the 'place hot' or cold buttons and click+drag the mouse around the picture to make more heat sources. There's lots of option! I went a little nuts.

The actual FEA code is very short, and is in the box below. Here's the description of how it works:

So the most important variables are 'now' and 'next', which are arrays filled with every temperature in our simulation. (just old and new from before).  At the start of our round they're identical. The outermost cells of these arrays are 'insulator', and will always have the same temperature as the cell it touches just inwards from them  (just like in real life).  For all of the cells representing Gallium (cells from 1 to numCells in both x and y) we need to calculate the temperature change each round. First we doublecheck the cell is gallium and not a heat source by checking the array isSource. If isSource[x][y] is false, we know the cell is gallium and we should continue, otherwise it's a heat source and we ignore it. For Gallium, we look in all 4 directions around our cell (up, down, left, right) and for each of these directions do the following:

change in temperature = coefficient * temperature difference.

When you add the changes from all 4 directions together, it becomes this line:

next[i][j] += coef*(now[i-1][j] + now[i+1][j] + now[i][j-1] + now[i][j+1] - 4*now[i][j]);

After doing this for all the cells, the gallium will be complete. Then, we go along the 4 edges of the 'next' array and set each insulator cell there to be the same temperature as the gallium cell it touches. Again, that's just how insulators works- no heat flows when there's no temperature difference! The index for insulator on the left side of the array is x=0, and the right side of the array is x=numCells+1. Top and bottom are likewise for y. To allow for these extra insulator cells, all of our arrays are actually numCells+2 large!

Now, the 'next' array has exactly what our temperatures should be for the next round. So, we step INTO THE FUTURE by copying 'next' back onto 'now', and increasing the system time by one tick.

The variables:

Spoiler

  // isSource: is a 2D arraw with an element for each cell. If a cell is a heat/cold source, isSource is "true"!
  // now: is a 2D arraw with an element for each cell. It contains the cell's current temperature as a float (a JS float)
  // next: is also a 2D arraw with an element for each cell. It contains the cell's temperature for the next time step
  // numCells: is the number of cells per row or collumn for our square domain. it's 40 in the link. The arrays are 2 larger than this to allow for insulator cells
  // coef: The transfer coefficient based on alpha

 

  for (i=1; i<=numCells; i++) {         //update temperatures
    for (j=1; j<=numCells; j++) {
      if (isSource[i][j] == false) {
	    next[i][j] += coef*(now[i-1][j] + now[i+1][j] + now[i][j-1] + now[i][j+1] - 4*now[i][j]);
      }
    }
  }

  for (i=1; i<=numCells; i++) {         //maintain temperature at insulator boundry conditions
    next[0][i] = next[1][i]; //left
    next[numCells+1][i] = next[numCells][i]; //right
    next[i][0] = next[i][1]; //top
    next[i][numCells+1] = next[i][numCells]; //bottom
  }
  
  for (i=0; i<=numCells+1; i++) {      //Copy 'next' onto 'now' to prepare for the next timestep. 
    for (j=0; j<=numCells+1; j++) {
      now[i][j] = next[i][j];
    }
  }
  
  curTime += timeStep;

 

 

Here's full the code. You can copy-paste it into a text editor, save it as an HTML file and run it from your own computer! Not sure why you would, but I need to have some excuse for using a language as knotty as JS. The 'important' FEA part is in a function called AdvanceFEA() { , and the active part of that is shown above in its entirety.

Spoiler

The commenting is terrible for the rest of this program- it's not really intended for other eyes. Just a bunch of UI junk anyways! I'm happy to answer questions, but just know in advance you'll find a lot of funny stuff in here if you look. Read at your own peril, you know?

<!DOCTYPE html>
<html lang="en">
<head>
  <meta charset="utf-8"/>
  <title>Simplest FEA</title>
</head>
<body>



<canvas id="displayCanvas" width="400" height="400" style="border:4px solid #3C3C3C;"></canvas>
<canvas id="bufferCanvas" width="100" height="100" style="border:4px solid #3C3C3C; display:none;"></canvas>

<h2 id="Monitor1"></h2>
<h2 id="Monitor2"></h2>

<button type="button" onclick='document.getElementById("Monitor1").innerHTML = "You clicked me! The HTML is still responding"'>Click Me!</button>
<button type="button" id="stepButton">Step</button>
<button type="button" id="runButton">Run</button>
<button type="button" id="stopButton" >Stop</button>
<button type="button" id="resetButton">Reset</button>

<br>
<button type="button" id="placeHotButton">Place Hot</button>
<button type="button" id="placeColdButton">Place Cold</button>
<button type="button" id="eraseButton">Erase</button>
<button type="button" id="clearAllButton">Clear All Sources</button>
<br><br>

<button type="button" id="changeTimeStepButton">Change Time Step</button>
Time Step (s): <textarea id="timeStepInput" rows="1" cols="10">0.002</textarea>
FEA steps/s = 60Hz * <textarea id="stepsPerFrameInput" rows="1" cols="10">59</textarea>
Redraws/s = 60Hz / ( 1 + <textarea id="redrawSkipsInput" rows="1" cols="10">2</textarea> )
<br><br>

<button type="button" id="reinitializeDomainButton">Reinitialize Domain</button>
Domain Size (in m): <textarea id="domainSizeInput" rows="1" cols="10">0.20</textarea>
Pixels Per Cell: <textarea id="pixelsPerCellInput" rows="1" cols="10">10</textarea>
Cells per row/col: <textarea id="numCellsInput" rows="1" cols="10">40</textarea>
<br><br>

<button type="button" id="setTempsButton">Set Temps</button>
Hot Temp (C): <textarea id="hotTempInput" rows="1" cols="10">50</textarea>
Cold Temp (C): <textarea id="coldTempInput" rows="1" cols="10">0</textarea>
Initial Temp (C): <textarea id="initialTempInput" rows="1" cols="10">25</textarea>
<br><br>

<button type="button" id="changeAlphaButton">Change Alpha</button>
Alpha (Thermal Diffusivity in m^2/s): <textarea id="alphaInput" rows="1" cols="10">0.00001855</textarea>

<h2 id="Monitor3"></h2>

<script>

"use strict";
//run rate, frame rate

var ctx = document.getElementById("displayCanvas").getContext("2d");
var buf = document.getElementById("bufferCanvas").getContext("2d");
var mainLoopRequest; //type: Object 

var runState = "stop";
var stepsPerFrame = 59; // multiply by 60 to get steps/second
var timeStep = 0.002; // in s
var redrawSkips = 2;

var curRedrawSkips = 0;
var curSteps = 0;

var domainSize = 0.20;  // in m
var pixelsPerCell_global = 10;
var numCells_global = 40;
var alpha_global = 0.00001855;
var coefficient_global = alpha_global*timeStep/((domainSize/numCells_global)*(domainSize/numCells_global));

var hotTemp_global = 50.00;
var coldTemp_global = 0.00;
var initialTemp = 25.00;

var img_global = buf.createImageData(numCells_global,numCells_global);
ctx.canvas.width = numCells_global*pixelsPerCell_global;
ctx.canvas.height = numCells_global*pixelsPerCell_global;
ctx.setTransform(pixelsPerCell_global, 0, 0, pixelsPerCell_global, 0, 0);
ctx.fillStyle = 'rgb(255, 165, 0)';
ctx.fillRect(0, 0, 150, 75);

ctx.imageSmoothingEnables = false;
buf.imageSmoothingEnables = false;

var placeType = "hot";
var mouseIsDown = "false";
var mouseOverCellX = 0;
var mouseOverCellY = 0;

var curTime = 0.00;

var i = 0;
var j = 0;
var k = 0;
var p = 0.00;
var q = 0.00;
var r = 0.00;

var F = 0;
var sumFEALoad = 0;
var feaLoadMonitor = [];
for (i=0; i<180; i++) feaLoadMonitor[i] = 0;

var D = 0;
var sumDisplayLoad = 0;
var displayLoadMonitor = [];
for (i=0; i<180; i++) displayLoadMonitor[i] = 0;


var now_global = [];
for (i=0; i<numCells_global+2; i++) {
  now_global[i] = [];
  for (j=0; j<numCells_global+2; j++) {
    now_global[i][j] = initialTemp;
  }
}

var next_global = [];
for (i=0; i<numCells_global+2; i++) {
  next_global[i] = [];
  for (j=0; j<numCells_global+2; j++) {
    next_global[i][j] = initialTemp;
  }
}

var isSource_global = [];
for (i=0; i<numCells_global+2; i++) {
  isSource_global[i] = [];
  for (j=0; j<numCells_global+2; j++) {
    isSource_global[i][j] = false;
  }
}


for (i=0; i<=10; i++) { //make demo heat sources
  now_global[15+i][1]  = hotTemp_global;
  next_global[15+i][1] = hotTemp_global;
  isSource_global[15+i][1] = true;

  now_global[25+i][15]  = hotTemp_global;
  next_global[15+i][15] = hotTemp_global;
  isSource_global[15+i][15] = true;  

  now_global[15+i][25]  = hotTemp_global;
  next_global[15+i][25] = hotTemp_global;
  isSource_global[15+i][25] = true;
  
  now_global[15+i][40]  = hotTemp_global;
  next_global[15+i][40] = hotTemp_global;
  isSource_global[15+i][40] = true;
  
  now_global[1][15+i]  = coldTemp_global;
  next_global[1][15+i] = coldTemp_global;
  isSource_global[1][15+i] = true;
  
  now_global[40][15+i]  = coldTemp_global;
  next_global[40][15+i] = coldTemp_global;
  isSource_global[40][15+i] = true;
  
  
}


/*
        isSource_global[x][y] = true;
		now_global[x][y]  = hotTemp_global;
		next_global[x][y] = hotTemp_global;
*/


function AdvanceFEA() {
  var isSource = isSource_global;
  var now = now_global; 
  var next = next_global;
  var numCells = numCells_global;
  var coef = coefficient_global;

  var i = 0;
  var j = 0;
  
  for (i=1; i<=numCells; i++) {
    for (j=1; j<=numCells; j++) {
      if (isSource[i][j] == false) {
	    next[i][j] += coef*(now[i-1][j] + now[i+1][j] + now[i][j-1] + now[i][j+1] - 4*now[i][j]);
	  }
    }
  }
  
  for (i=1; i<=numCells; i++) { //maintain temperature at insulator boundry conditions
    next[0][i] = next[1][i];
    next[numCells+1][i] = next[numCells][i];
    next[i][0] = next[i][1];
    next[i][numCells+1] = next[i][numCells];
  }
  
  for (i=0; i<=numCells+1; i++) { //
    for (j=0; j<=numCells+1; j++) {
      now[i][j] = next[i][j];
    }
  }
  
  curTime += timeStep;
}



//Redraw(); //Run once on startup NO LONGER
function Redraw() { 
  
  //make locals (copies for scalars and references for arrays)
  var numCells = numCells_global; 
  var pixelsPerCell = pixelsPerCell_global;
  var middleTemp = (hotTemp_global + coldTemp_global)/2;
  var scale = 255/(hotTemp_global - coldTemp_global);
  var now = now_global; //passes an array reference for the Temperatures
  var id = img_global.data; //reference to a clamped Uint8 array for display buffering
//  var id = uintc8;
 
 
  var curX = 0;
  var curY = 0;
  
  var r = 0; //red
  var g = 0; //green
  var b = 0; //blue
  
  var i = 0;
  var j = 0;
  var k = 0;

  k=0;
  for (i=1; i<=numCells; i++) { //skip the outermost cells: 0 and numCells+1
    for (j=1; j<=numCells; j++) {
	  //id auto clamps and rounds
	  id[k++] = 127.5+scale*(now[j][i]-middleTemp);   //red
	  id[k++] = 127.5-scale*Math.abs(now[j][i]-middleTemp);  //green
	  id[k++] = 127.5-scale*(now[j][i]-middleTemp);  //blue
	  id[k++] = 255;
    // the k++s read the current index then increment, totaling k+=4;
    } 
  }
  
 //   console.log(id);
  
  buf.putImageData(img_global,0,0);
  ctx.drawImage(buf.canvas,0,0); //ctx is autoscaled, but otherwise identical to buf.

  document.getElementById("Monitor1").innerHTML = "Cell (" + mouseOverCellX + "," + mouseOverCellY + ") Temp = " + now_global[mouseOverCellX][mouseOverCellY].toFixed(12) + " C"; 
  
  var feaLoad = sumFEALoad/(.01*180*16.667); //sumFEALoad is a walking average in essence, it's found in the MainLoop function. We just need to scale it into % now.
  var displayLoad = sumDisplayLoad/(.01*180*16.667);

  document.getElementById("Monitor2").innerHTML = "Current Time = " + curTime.toFixed(5) + " s <br>  FEA Load = " + feaLoad.toFixed(0) + "%    Display Load = " + displayLoad.toFixed(0) + "%";
  
  document.getElementById("Monitor3").innerHTML = "Coefficient = " + coefficient_global.toFixed(10) + " (unitless) ";
  
  //console.log(displayTime);
}

requestAnimationFrame(MainLoop); //run once on startup, noncancellable
function MainLoop() {

  var time1 = performance.now();

  var i=0;  
  for (i=0; i<stepsPerFrame; i++) AdvanceFEA();
  
  
  var time2 = performance.now();  
  
  curRedrawSkips++; //global scope
  if (curRedrawSkips >= redrawSkips) {
    curRedrawSkips = 0;
    Redraw();
  }
  //the next line runs this once per frame. Very handy for under the hood reasons. mainLoopRequest provides a way to cancel the request.
  if (runState = "run") mainLoopRequest = requestAnimationFrame(MainLoop);

  var time3 = performance.now();
  
  
  //Load Monitoring
  var feaTime = (time2 - time1);
  sumFEALoad += feaTime - feaLoadMonitor[F]; //this is a walking average, so pick up a new value and drop the old one.
  feaLoadMonitor[F] = feaTime;
  F++;
  if (F>=180) F=0;
  
  var displayTime = (time3 - time2);
  sumDisplayLoad += displayTime - displayLoadMonitor[D]; //this is a walking average, so pick up a new value and drop the old one.
  displayLoadMonitor[D] = displayTime;
  D++;
  if (D>=180) D=0;  
  
}




function ChangeAlpha() {
  alpha_global = parseFloat(document.getElementById("alphaInput").value);
  coefficient_global = alpha_global*timeStep/((domainSize/numCells_global)*(domainSize/numCells_global));
  Redraw();
  
}


function ReinitializeDomain() {
   
  runState = "stop";
  domainSize = parseFloat(document.getElementById("domainSizeInput").value); 
  pixelsPerCell_global = parseInt(document.getElementById("pixelsPerCellInput").value);
  numCells_global = parseInt(document.getElementById("numCellsInput").value);
  coefficient_global = alpha_global*timeStep/((domainSize/numCells_global)*(domainSize/numCells_global));
  
  
  buf.canvas.width = numCells_global;
  buf.canvas.height = numCells_global;
  img_global = buf.createImageData(numCells_global,numCells_global);
  ctx.canvas.width = numCells_global*pixelsPerCell_global;
  ctx.canvas.height = numCells_global*pixelsPerCell_global;
  ctx.setTransform(pixelsPerCell_global, 0, 0, pixelsPerCell_global, 0, 0);
  
  mouseOverCellX = Math.ceil(numCells_global/2);
  mouseOverCellY = Math.ceil(numCells_global/2);


  var i = 0; //locals
  var j = 0;
  
  now_global = [];
  for (i=0; i<numCells_global+2; i++) {
    now_global[i] = [];
    for (j=0; j<numCells_global+2; j++) {
      now_global[i][j] = initialTemp;
    }
  }

  next_global = [];
  for (i=0; i<numCells_global+2; i++) {
    next_global[i] = [];
    for (j=0; j<numCells_global+2; j++) {
      next_global[i][j] = initialTemp;
    }
  }  
  
  isSource_global = [];
  for (i=0; i<numCells_global+2; i++) {
    isSource_global[i] = [];
    for (j=0; j<numCells_global+2; j++) {
      isSource_global[i][j] = false;
    }  
  }

  
  curTime = 0.00;
  Redraw();
  
}

function ChangeTimeStep() {

  timeStep = parseFloat(document.getElementById("timeStepInput").value);
  stepsPerFrame = parseInt(document.getElementById("stepsPerFrameInput").value);
  redrawSkips = parseInt(document.getElementById("redrawSkipsInput").value);

  coefficient_global = alpha_global*timeStep/((domainSize/numCells_global)*(domainSize/numCells_global));
  
  Redraw();
  
}

function SetTemps() {

  hotTemp_global = parseFloat(document.getElementById("hotTempInput").value);
  coldTemp_global = parseFloat(document.getElementById("coldTempInput").value);
  initialTemp = parseFloat(document.getElementById("initialTempInput").value);
  
  Redraw();
  
}

function HandleMouseEvents(e) {
  mouseOverCellX = Math.floor(e.offsetX/pixelsPerCell_global)+1; //+1 because cells display starting from 1 through numCells
  mouseOverCellY = Math.floor(e.offsetY/pixelsPerCell_global)+1;
  var x = mouseOverCellX;
  var y = mouseOverCellY;
  

  if (x>=1 && y>=1 && x<=numCells_global && y<=numCells_global) {    //if in bounds...
    document.getElementById("Monitor1").innerHTML = "Cell (" + x + "," + y + ") Temp = " + now_global[x][y].toFixed(12) + " C";
	
    if (mouseIsDown == true) {
      if (placeType == "erase") {
        isSource_global[x][y] = false;
		now_global[x][y]  = initialTemp;
		next_global[x][y] = initialTemp;
      }	  
	  else if (placeType == "hot") {
        isSource_global[x][y] = true;
		now_global[x][y]  = hotTemp_global;
		next_global[x][y] = hotTemp_global;
	  }
	  else if (placeType == "cold") {
        isSource_global[x][y] = true;
		now_global[x][y]  = coldTemp_global;
		next_global[x][y] = coldTemp_global;
	  }
	  Redraw();
    }		
	
  }
}

//Button Listeners
//Run
document.getElementById("stepButton").onclick = function() {
  runState = "stop";
  cancelAnimationFrame(mainLoopRequest);
  AdvanceFEA();
  Redraw();
};

document.getElementById("runButton").onclick = function() {
  if (runState != "run") {
    runState = "run";
    mainLoopRequest = requestAnimationFrame(MainLoop); //runs MainLoop
  }
};

//Stop
document.getElementById("stopButton").onclick = function() {
  runState = "stop";
  cancelAnimationFrame(mainLoopRequest)
};

//Reset
document.getElementById("resetButton").onclick = function() {
  
  var isSource = isSource_global;
  var now = now_global; 
  var next = next_global;
  var numCells = numCells_global;

  var i = 0;
  var j = 0;
  for (i=1; i<=numCells; i++) {
    for (j=1; j<=numCells; j++) {
      if (isSource[i][j] == false) {
	    now[i][j] = initialTemp;
	    next[i][j] = initialTemp;
	  }
    }
  }
  
  curTime = 0.00;
  Redraw();
 
};

//Clear All
document.getElementById("clearAllButton").onclick = function() {
  var isSource = isSource_global;
  var now = now_global; 
  var next = next_global;
  var numCells = numCells_global;

  var i = 0;
  var j = 0;
  for (i=1; i<=numCells; i++) {
    for (j=1; j<=numCells; j++) {
      isSource[i][j] = false;
    }
  }
  Redraw();
  
};

//Place Hot
document.getElementById("placeHotButton").onclick = function() {
  placeType = "hot";
};

//Place Cold
document.getElementById("placeColdButton").onclick = function() {
  placeType = "cold";
};

//Erase
document.getElementById("eraseButton").onclick = function() {
  placeType = "erase";
};

//Change Timesteps
document.getElementById("changeTimeStepButton").onclick = function() {
  ChangeTimeStep();
};

//Reinitialize Domain
document.getElementById("reinitializeDomainButton").onclick = function() {
  ReinitializeDomain();
};

//Set Temps
document.getElementById("setTempsButton").onclick = function() {
  SetTemps();
};

document.getElementById("changeAlphaButton").onclick = function() {
  ChangeAlpha();
};



//Mouse Listeners
document.getElementById('displayCanvas').onmousedown = function(e) {
  mouseIsDown = true;
  HandleMouseEvents(e); //e contains details on the mouse event
};
document.onmouseup = function() {
  mouseIsDown = false;
};

document.getElementById('displayCanvas').onmousemove = function(e) {
  HandleMouseEvents(e); //e contains details on the mouse event
};

document.getElementById('displayCanvas').onmouseleave = function(e) {
  mouseOverCellX = Math.ceil(numCells_global/2);
  mouseOverCellY = Math.ceil(numCells_global/2); 
  document.getElementById("Monitor1").innerHTML = "Cell (" + mouseOverCellX + "," + mouseOverCellY + ") Temp = " + now_global[mouseOverCellX][mouseOverCellY].toFixed(12) + " C";  
};


//document.getElementById('displayCanvas').onmouseenter = divOnMouseenter;


</script>



</body>
</html>

 

Alright, time to disappear for a while and pass this torch along. Cheers!

Edited June 12, 2019 by Cunjo Carl