% stefan2D.prb P t1 4000000.0 P dt 20000.0 P L 334000 % latent heat (ice, mks units) P k 2.26 % thermal diffusivity P rho 1000. % density P cp 4200. % heat capacity P T0 0.01 A heat e V400 = [grad]{V200}_j e V600 = minv V400 e V800 = V601,V603,V602,V604 e V1000 = V800 P V600 e V1201 = (V401*V404-V402*V403) e V1300 = rho*cp*(V301*V1201,V302*V1201) e {ondt*V1201}U1-{V600 p V1300}_jD_jU1 = \ {ondt*V1201}{V101}+D_j{k*V1001*V1201, \ k*V1002*V1201,k*V1003*V1201, \ k*V1004*V1201}_jk D_k U1 b 1 U1 = {V101} b 2 U1 = {V101} A spread e V400 = [grad] {V200}_j e V600 = minv V400 e V800 = V601,V603,V602,V604 e V1000 = V800 P V600 e V300_i = [grad] {V101}_i e D_j D_j U1_i + D_j D_i U1_j b 2 {L*rho*V600}_jU1_j-{V1000}_j{k*V300}_j={0.0} b 1 U1_1 = 0. b 3 U1_2 = 0. < run t=0 ondt=rho*cp/dt V101=(T0+X2*X2)*(X1-1) % temperature V200=X1*.02,X2 % mapping fn x=x(X,t) V300=0.,0. % time deriv dx/dt while t<t1 heat 1 100 V100 200 V200 300 V300 solve V201=V101 popp spread solve V400=V100 V300=V300+V400*dt popp mapm 200 black contour 101 mapm 200 t=t+dt endwhile >