Minimization of Free Energy

In the section Description of Introductory Example the description of the steady-state heat conduction on a three-dimensional domain was given. The solution of the same physical problem can be obtained also as a minimum of the free energy of the problem. Free energy of the heat conduction problem can be formulated as
Π=(∫_Ω)((1/2)k Δϕ.Δϕ-ϕ Q)dΩ  
where a ϕ indicates temperature, a k is the conductivity and a Q is the heat generation per unit volume and Ω is the domain of the problem.

The domain of the example is a cube filled with water ([-.0.5m,0.5m]×[-0.5m,0.5m]×[0,1m]). On all sides, apart from the upper surface, the constant temperature ϕ=0 is maintained. The upper surface is isolated so that there is no heat flow over the boundary. There exists a constant heat source Q=500 W/(m^3) inside the cube. The thermal conductivity of water is 0.58 W/m K. The task is to calculate the temperature distribution inside the cube.

The problem is formulated using various approaches:
A.  Trial polynomial interpolation
	M.G	Gradient method of optimization + Mathematica directly
	M.N	Newton method of optimization + Mathematica directly
	A.G	Gradient method of optimization + AceGen+MathLink
	A.N	Newton method of optimization + AceGen+MathLink
B.	Finite difference interpolation 
	M.G	Gradient method of optimization + Mathematica directly
	M.N	Newton method of optimization + Mathematica directly
	A.G	Gradient method of optimization + AceGen+MathLink
	A.N	Newton method of optimization + AceGen+MathLink
C.AceFEM	Finite element method

The following quantities are compared:
•	temperature at the central point of the cube (ϕ(0.,0.,0.5))
•	time for derivation of the equations
•	time for solution of the optimization problem
•	number of unknown parameters used to discretize the problem
•	peak memory allocated during the analysis
•	number of evaluations of function, gradient and hessian.

({{Method, mesh,  ϕ, derivation
time(s), solution
time (s), No. of variables, memory (MB), No. of calls}, {A.MMA.Gradient , 5×5×5  , 55.9, 8.6, 56.0, 80, 136, 964}, {A.MMA.Newton, 5×5×5  , 55.9, 8.6, 2588.3, 80, 1050, 4}, {A.AceGen.Gradient, 5×5×5  , 55.9, 6.8, 3.3, 80, 4, 962}, {A.AceGen.Newton, 5×5×5  , 55.9, 13.0, 0.8, 80, 4, 4}, {B.MMA.Gradient, 11×11×11  , 57.5, 0.3, 387.5, 810, 10, 1685}, {B.MMA.Newton, 11×11×11 , 57.5, 0.3, 4.2, 810, 16, 4}, {B.AceGen.Gradient, 11×11×11 , 57.5, 1.4, 28.16, 810, 4, 1598}, {B.AceGen.Newton, 11×11×11 , 57.5, 4.0, 1.98, 810, 4, 4}, {C.AceFEM, 10×10×10, 56.5, 5.0, 2.0, 810, 6, 2}, {C.AceFEM, 20×20×20, 55.9, 5.0, 3.2, 7220, 32, 2}, {C.AceFEM, 30×30×30, 55.9, 5.0, 16.8, 25230, 139, 2}})

The case A with the trial polynomial interpolation represents the situation where the merit function is complicated and the number of parameters is small. The case B with the  finite difference interpolation represents the situation where the merit function is simple and the number of parameters is large. 
REMMARK: The presented example is meant to illustrate the general symbolic approach to minimization of complicated merit functions and is not the state of the art solution of thermal conduction problem.

A) Trial Lagrange polynomial interpolation

A trial function for temperature ϕ is constructed as a fifth order Lagrange polynomial in x y and z direction. The chosen trial function is constructed in a way that  satisfies boundary conditions.

<<AceGen`;
Clear[x,y,z,α];
kcond=0.58;Q=500;
order=5;
nterm=(order-1)(order-1)(order)

Here the fifth order Lagrange polynomials are constructed in three dimensions.

toc=Table[{x,0},{x,-0.5,0.5,1/order}];xp=MapIndexed[ InterpolatingPolynomial[ReplacePart[toc,1,{#,2}],x]&,Range[2,order]];
yp=MapIndexed[ InterpolatingPolynomial[ReplacePart[toc,1,{#,2}],y]&,Range[2,order]];
toc=Table[{x,0},{x,0.,1.,1/order}];
zp=MapIndexed[ InterpolatingPolynomial[ReplacePart[toc,1,{#,2}],z]&,Range[2,order+1]];
ϕi=Array[α,nterm];
poly=Flatten[Outer[Times,xp,yp,zp]]//Chop;
ϕ=poly.ϕi;

poly[[28]]
Plot3D[poly[[28]]/.z->0.5,{x,-0.5,0.5},{y,-0.5,0.5},PlotRange->All]

Here the Gauss points and weights are calculated for ngp×ngp×ngp Gauss numerical integration of the free energy over the domain [-.0.5m,0.5m]×[-0.5m,0.5m]×[0,1m].

ngp=6;
<<NumericalDifferentialEquationAnalysis`;
g1=GaussianQuadratureWeights[ngp,-0.5,0.5];
g2=GaussianQuadratureWeights[ngp,-0.5,0.5];
g3=GaussianQuadratureWeights[ngp,0,1];
gp={g1[[#1[[1]],1]],g2[[#1[[2]],1]],g3[[#1[[3]],1]],g1[[#1[[1]],2]]*g2[[#1[[2]],2]]*g3[[#1[[3]],2]]}&/@Flatten[Array[{#3,#2,#1}&,{ngp,ngp,ngp}],2];

The subsection Definitions has to be executed before the current subsection.

start=SessionTime[];
Δϕ={D[ϕ,x],D[ϕ,y],D[ϕ,z]};
Π=1/2 kcond Δϕ.Δϕ-ϕ Q;
Πi=Total[Map[(#[[4]] Π/.{x->#[[1]],y->#[[2]],z->#[[3]]})&,gp]];
derivation=SessionTime[]-start

G. Gradient based optimization

start=SessionTime[];ii=0;
sol=FindMinimum[Πi,Array[{α[#] ,0.}&,nterm],Method->"Gradient",EvaluationMonitor:>(ii++;)];
{ii,ϕ/.sol[[2]]/.{x->0,y->0,z->0.5}}
SessionTime[]-start

N. Newton method based optimization

start=SessionTime[];ii=0;
sol=FindMinimum[Πi,Array[{α[#] ,0.}&,nterm],Method->"Newton",EvaluationMonitor:>(ii++;)];
{ii,ϕ/.sol[[2]]/.{x->0,y->0,z->0.5}}
SessionTime[]-start

The subsection Definitions has to be executed before the current subsection.

start=SessionTime[];SMSInitialize["Thermal","Environment"->"MathLink","Mode"->"Prototype","ADMethod"->"Forward"]

Πf[i_]:=(
ai⊨SMSReal[Array[a$$,nterm]];
ag⊨SMSArray[ai];
{xa,ya,za,wa}⊨Map[SMSArray,Transpose[gp]];
{xi,yi,zi}⊨SMSFreeze[{SMSPart[xa,i],SMSPart[ya,i],SMSPart[za,i]}];
{xpr,ypr,zpr}⊨{xp/.x->xi,yp/.y->yi,zp/.z->zi};
poly⊨SMSArray[Flatten[Outer[Times,xpr,ypr,zpr]]];
ϕt⊨SMSDot[poly,ag];
Δϕ⊨SMSD[ϕt,{xi,yi,zi}];
wi⊨SMSPart[wa,i];
wi(1/2 kcond Δϕ.Δϕ-ϕt Q)
)

SMSModule["FThermal",Real[a$$[nterm],f$$],"Input"->a$$,"Output"->f$$];
SMSExport[0,f$$];
SMSDo[i,1,gp//Length];
  Π⊨Πf[i];
  SMSExport[ Π,f$$,"AddIn"->True];
SMSEndDo[];

SMSModule["GThermal",Real[a$$[nterm],g$$[nterm]],"Input"->a$$,"Output"->g$$];
SMSExport[Table[0,{nterm}],g$$];
SMSDo[i,1,gp//Length];
  Π⊨Πf[i];
  SMSDo[j,1,nterm];
    δΠ⊨SMSD[Π,ag,j,"Method"->"Forward"];
    SMSExport[δΠ,g$$[j],"AddIn"->True];
  SMSEndDo[];
SMSEndDo[];

derivation=SessionTime[]-start

SMSModule["HThermal",Real[a$$[nterm],h$$[nterm,nterm]],"Input"->a$$,"Output"->h$$];
SMSDo[i,1,nterm];
  SMSDo[j,1,nterm];
    SMSExport[0,h$$[i,j]];
  SMSEndDo[];
SMSEndDo[];
SMSDo[i,1,gp//Length];
  Π⊨Πf[i];
  SMSDo[j,1,nterm];
    δΠ⊨SMSD[Π,ag,j,"Method"->"Forward"];
    SMSDo[k,1,nterm];
      hij⊨SMSD[δΠ,ag,k,"Method"->"Forward"];
      SMSExport[hij,h$$[j,k],"AddIn"->True];
    SMSEndDo[];
  SMSEndDo[];
SMSEndDo[];

SMSWrite[];

SMSInstallMathLink["Optimize"->False]
derivation=SessionTime[]-start

G. Gradient based optimization

start=SessionTime[];ii=0;
sol=FindMinimum[FThermal[ϕi],{ϕi,Table[0,{nterm}]},Method->"Gradient",Gradient->GThermal[ϕi],EvaluationMonitor:>(ii++;)];
{ii,ϕ/.MapThread[Rule,List@@sol[[2,1]]]/.{x->0,y->0,z->0.5},SessionTime[]-start}

N. Newton method based optimization

start=SessionTime[];ii=0;
sol=FindMinimum[FThermal[ϕi],{ϕi,Table[0,{nterm}]},
Method->{"Newton", Hessian->HThermal[ϕi]},
Gradient->GThermal[ϕi],EvaluationMonitor:>(ii++;)];
{ii,ϕ/.MapThread[Rule,List@@sol[[2,1]]]/.{x->0,y->0,z->0.5},SessionTime[]-start}

B) Finite difference interpolation

The central difference approximation of derivatives is used for the points inside the cube and backward or forward difference for the points on the boundary.

<<AceGen`;
Clear[α,i,j,k];
nx=ny=nz=11;
dlx=1./(nx-1);
dly=1./(ny-1);
dlz=1./(nz-1);
bound={0};
nboun=1;
kcond=0.58;Q=500;

nterm=0;dofs={};
index=Table[Which[
i≤2 || i ≥  nx+1  || j≤2 || j ≥  ny+1  || k≤ 2 ,b[1]
,k==nz+2,
If[FreeQ[dofs,α[i,j,k-1]]
,++nterm;AppendTo[dofs,α[i,j,k-1]->nterm];nterm
,α[i,j,k-1]/.dofs
]
,True,
If[FreeQ[dofs,α[i,j,k]]
,++nterm;AppendTo[dofs,α[i,j,k]->nterm];nterm
,α[i,j,k]/.dofs
]
],
{i,1,nx+2},{j,1,ny+2},{k,1,nz+2}]/.b[i_]:>nterm+i;
ϕi=Array[α,nterm];
nterm

The subsection Definitions have to be executed before the current subsection.

start=SessionTime[];
Πi=Sum[
dlxt=If[i==2 || i==nx+1 ,dlxt=dlx/2,dlx];
dlyt=If[j==2 || j==ny+1 ,dlyt=dly/2,dly];
dlzt=If[k==2 || k==nz+1 ,dlzt=dlz/2,dlz];
vol=dlxt dlyt dlzt;
aijk=Map[If[#>nterm,bound[[#-nterm]],α[#]]&,
Extract[index,{{i,j,k},{i-1,j,k},{i+1,j,k},{i,j-1,k},
{i,j+1,k},{i,j,k-1},{i,j,k+1}}]];
grad={((aijk[[3]]-aijk[[2]])/(2 dlxt)),((aijk[[5]]-aijk[[4]])/(2 dlyt)),((aijk[[7]]-aijk[[6]])/(2 dlzt))};

vol(1/2 kcond grad.grad-Q aijk[[1]])
,{i,2,nx+1},{j,2,ny+1},{k,2,nz+1}
];
derivation=SessionTime[]-start

G. Gradient based optimization

start=SessionTime[];ii=0;
sol=FindMinimum[Πi,Array[{α[#] ,0.}&,nterm],Method->"Gradient",EvaluationMonitor:>(ii++;)];
{ii,α[index[[(nx+3)/2,(ny+3)/2,(nz+3)/2]]]/.sol[[2]],SessionTime[]-start}

N. Newton method based optimization

start=SessionTime[];ii=0;
sol=FindMinimum[Πi,Array[{α[#] ,0.}&,nterm],Method->"Newton",EvaluationMonitor:>(ii++;)];
{ii,α[index[[(nx+3)/2,(ny+3)/2,(nz+3)/2]]]/.sol[[2]],SessionTime[]-start}

The subsection Definitions have to be executed before the current subsection.

start=SessionTime[];SMSInitialize["Thermal","Environment"->"MathLink","Mode"->"Prototype","ADMethod"->"Backward"]

Πf[i_,j_,k_]:=(
indexp⊨SMSInteger[Map[
index$$[(#[[1]]-1)*(nyp+2)(nzp+2)+(#[[2]]-1)*(nzp+2)+#[[3]]]&,
{{i,j,k},{i-1,j,k},{i+1,j,k},{i,j-1,k},
{i,j+1,k},{i,j,k-1},{i,j,k+1}}]];
aijk⊨SMSReal[Map[a$$[#]&,indexp]];
{dx,dy,dz,kc,Qt}⊨SMSReal[Array[mc$$,5]];
SMSIf[i==2 || i==nxp+1 ];
  dlxt⫤dx/2;
SMSElse[];
  dlxt⊣dx;
SMSEndIf[dlxt];
SMSIf[j==2 || j==nyp+1 ];
  dlyt⫤dy/2;
SMSElse[];
  dlyt⊣dy;
SMSEndIf[dlyt];
SMSIf[k==2 || k==nzp+1 ];
  dlzt⫤dz/2;
SMSElse[];
  dlzt⊣dz;
SMSEndIf[dlzt];
vol⊨dlxt dlyt dlzt;
grad={((aijk[[3]]-aijk[[2]])/(2 dlxt)),((aijk[[5]]-aijk[[4]])/(2 dlyt)),((aijk[[7]]-aijk[[6]])/(2 dlzt))};
vol (1/2 kc grad.grad-Qt aijk[[1]])
)

SMSModule["FThermal",Integer[ndof$$,nt$$[3],index$$["*"]],Real[a$$["*"],mc$$["*"],f$$],"Input"->{ndof$$,nt$$,index$$,a$$,mc$$},"Output"->f$$];
SMSExport[0,f$$];
{nxp,nyp,nzp}⊨SMSInteger[Array[nt$$,3]];
SMSDo[i,2,nxp+1];
  SMSDo[j,2,nyp+1];
    SMSDo[k,2,nzp+1];
      Π⊨Πf[i,j,k];
      SMSExport[Π,f$$,"AddIn"->True];
    SMSEndDo[];
  SMSEndDo[];
SMSEndDo[];

SMSModule["GThermal",Integer[ndof$$,nt$$[3],index$$["*"]],Real[a$$["*"],mc$$["*"],g$$[ndof$$]],"Input"->{ndof$$,nt$$,index$$,a$$,mc$$},"Output"->g$$];
ndof⊨SMSInteger[ndof$$];
{nxp,nyp,nzp}⊨SMSInteger[Array[nt$$,3]];
SMSDo[i,1,ndof];
  SMSExport[0,g$$[i]];
SMSEndDo[];

SMSDo[i,2,nxp+1];
  SMSDo[j,2,nyp+1];
    SMSDo[k,2,nzp+1];
      Π⊨Πf[i,j,k];
      SMSDo[i1,1,indexp//Length];
        dof⊨SMSPart[indexp,i1];
        SMSIf[dof<=ndof];
          gi⊨SMSD[Π,aijk,i1];
          SMSExport[gi,g$$[dof],"AddIn"->True];
        SMSEndIf[];
      SMSEndDo[];
    SMSEndDo[];
  SMSEndDo[];
SMSEndDo[];

derivation=SessionTime[]-start

SMSModule["HThermal",Integer[ndof$$,nt$$[3],index$$["*"]],Real[a$$["*"],mc$$["*"],h$$[ndof$$,ndof$$]],"Input"->{ndof$$,nt$$,index$$,a$$,mc$$},"Output"->h$$];
ndof⊨SMSInteger[ndof$$];
{nxp,nyp,nzp}⊨SMSInteger[Array[nt$$,3]];
SMSDo[i,1,ndof];
  SMSDo[j,1,ndof];
    SMSExport[0,h$$[i,j]];
  SMSEndDo[];
SMSEndDo[];

SMSDo[i,2,nxp+1];
  SMSDo[j,2,nyp+1];
    SMSDo[k,2,nzp+1];
      Π⊨Πf[i,j,k];
      SMSDo[i1,1,indexp//Length];
        dofi⊨SMSPart[indexp,i1];
        SMSIf[dofi<=ndof];
          gi⊨SMSD[Π,aijk,i1];
          SMSDo[j1,1,indexp//Length];
            dofj⊨SMSPart[indexp,j1];
            SMSIf[dofj<=ndof];
              hij⊨SMSD[gi,aijk,j1];
              SMSExport[hij,h$$[dofi,dofj],"AddIn"->True];
            SMSEndIf[];
          SMSEndDo[];
        SMSEndIf[];
      SMSEndDo[];
    SMSEndDo[];
  SMSEndDo[];
SMSEndDo[];

SMSWrite[];

SMSInstallMathLink["Optimize"->True]
derivation=SessionTime[]-start

G. Gradient based optimization

start=SessionTime[];ii=0;
indexb=Flatten[index];
sol=FindMinimum[
FThermal[nterm,{nx,ny,nz},indexb,Join[ϕi,bound],{dlx,dly,dlz,kcond,Q}]
,{ϕi,Table[0,{nterm}]},
Method->"Gradient",
Gradient->GThermal[nterm,{nx,ny,nz},indexb,Join[ϕi,bound],{dlx,dly,dlz,kcond,Q}]
,EvaluationMonitor:>(ii++;)];{ii,α[index[[(nx+3)/2,(ny+3)/2,(nz+3)/2]]]/.MapThread[Rule,List@@sol[[2,1]]],SessionTime[]-start}

N. Newton method based optimization

start=SessionTime[];ii=0;
indexb=Flatten[index/.b[i_]:>nterm+i];
sol=FindMinimum[
FThermal[nterm,{nx,ny,nz},indexb,Join[ϕi,bound],{dlx,dly,dlz,kcond,Q}]
,{ϕi,Table[0,{nterm}]},
Method->{"Newton", Hessian->HThermal[nterm,{nx,ny,nz},indexb,Join[ϕi,bound],{dlx,dly,dlz,kcond,Q}]},
Gradient->GThermal[nterm,{nx,ny,nz},indexb,Join[ϕi,bound],{dlx,dly,dlz,kcond,Q}]
,EvaluationMonitor:>(ii++;)];{ii,α[index[[(nx+3)/2,(ny+3)/2,(nz+3)/2]]]/.MapThread[Rule,List@@sol[[2,1]]],SessionTime[]-start}

The tangent matrix is in the case of finite difference approximation extremely sparse.

MatrixPlot[HThermal[nterm,{nx,ny,nz},indexb,Join[0 ϕi,bound],{dlx,dly,dlz,kcond,Q}]]

C) Finite element method

First the finite element mesh 30×30×30 is used to obtain convergence solution at the central point of the cube. The procedure to generate heat-conduction element that is used in this example is explained in AceGen manual section Standard FE Procedure .

<<AceFEM`;
start=SessionTime[];
SMTInputData[];
k=0.58;Q=500;
nn=10;
SMTAddDomain["cube","ExamplesHeatConduction",{k,0,0,Q}];
SMTAddEssentialBoundary[{"X"==-0.5|| "X"==0.5 || "Y"==-0.5 || "Y"==0.5 || "Z"==0.&,0}];
SMTMesh["cube","H1",{nn,nn,nn},{
{{{-0.5,-0.5,0},{0.5,-0.5,0}},{{-0.5,0.5,0},{0.5,0.5,0}}},{{{-0.5,-0.5,1},{0.5,-0.5,1}},{{-0.5,0.5,1},{0.5,0.5,1}}}
}];
SMTAnalysis[];

SMTNextStep[0,1];
SMTNewtonIteration[];

SMTPointValues[{0,0,0.5},SMTPost[1]]
SessionTime[]-start

SMTShowMesh["Mesh"->True,"Elements"->True,"Field"->SMTPost[1]]