Venkat Subramanian

386 Reputation

13 Badges

14 years, 205 days

MaplePrimes Activity


These are questions asked by


 

NULL

restart:

Digits:=15;

15

A:=Matrix(4,4,[[-1,2,0,0],[2,-1,2,0],[0,2,-1,2],[0,0,1,-1]],datatype=float[8],storage=sparse);

Matrix(4, 4, {(1, 1) = -1.0000000000000000, (1, 2) = 2.0000000000000000, (1, 3) = 0., (1, 4) = 0., (2, 1) = 2.0000000000000000, (2, 2) = -1.0000000000000000, (2, 3) = 2.0000000000000000, (2, 4) = 0., (3, 1) = 0., (3, 2) = 2.0000000000000000, (3, 3) = -1.0000000000000000, (3, 4) = 2.0000000000000000, (4, 1) = 0., (4, 2) = 0., (4, 3) = 1.0000000000000000, (4, 4) = -1.0000000000000000})

b:=Vector(4,[1,0.,0,0],datatype=float[8]):

 

sol:=LinearAlgebra:-LinearSolve(A,b,method=SparseDirectMKL);

Error, invalid input: LinearAlgebra:-LinearSolve expects value for keyword parameter method to be of type identical(none,SparseLU,SparseDirect,SparseDirectMKL,SparseIterative,LU,QR,solve,hybrid,Cholesky,subs,modular), but received SparseDirectMKL

NULL


 

Download buglinearsolve.mw

I was pleased to see the description of SparseDirectMKL, but it is not implemented properly yet. It is a step in the right direction, so please make this available in the near future

 

In a recent question/conversation, I had discussed integrating dsolve/numeric-based codes with NLPsolve at 
https://www.mapleprimes.com/questions/236494-Inconsistent-Behavior-With-Dsolvenumeric-And-NLPSolve-sqp

@acer was able to help resolve the issue by calling NLPSolve with higher optimality tolerance.

I am starting a new question/post to show the need for evalhf as the previous post takes too long to load.
Code is attached for the same.
 

restart:

currentdir();

"C:\Users\Venkat\OneDrive - The University of Texas at Austin\Documents\trial\Learnningsqpoptim"

Test code written by Dr. Venkat Subramanian at UT Austin, 05/31/2023 and has been updated multiple times. This code uses CVP approach (piecwise constant) to perform optimal control. NLPSolve combination with RK2 approach to integrate ODEs (constant step size) works well for this problem. But dsolve numeric based codes work (with optimality tolerance fix from acer), but cannot handle large values of nvar (optimization variables).

Digits:=15;

15

 

eqodes:=[diff(ca(t),t)=-(u+u^2/2)*1.0*ca(t),diff(cb(t),t)=1.0*u*ca(t)-0.1*cb(t)];

[diff(ca(t), t) = -1.0*(u+(1/2)*u^2)*ca(t), diff(cb(t), t) = 1.0*u*ca(t)-.1*cb(t)]

soln:=dsolve({op(eqodes),ca(0)=alpha,cb(0)=beta},type=numeric,'parameters'=[alpha,beta,u],savebinary=true):

 

 

Note that Vector or Array form can be used for RK2h to implement the procedure for any number of variables/ODEs. But the challenge will be when implicit methods are used to integrate ODEs/DAEs and running them in evalhf form (this can be done with Gauss elimination type linear solver, but this will be limited to small number of ODEs/DAEs say 200 or so).

RK2h:=proc(NN,u,dt,y0::Vector(datatype=float[8]))#this procedure can be made efficient in vector form
local j,c1mid,c2mid,c10,c20,c1,c2;
option hfloat;
c10:=y0[1];c20:=y0[2];
for j from 1 to NN do
  c1mid:=c10+dt/NN*(-(u+u^2/2)*c10):
  c2mid:=c20+dt/NN*(u*c10)-0.1*dt/NN*c20:
  c1:=c10/2+c1mid/2+dt/NN/2*(-(u+u^2/2)*c1mid):
  c2:=c20/2+c2mid/2+dt/NN/2*(u*c1mid)-0.1*dt/NN/2*c2mid:
  c10:=c1:c20:=c2:
  od:
y0[1]:=c10:y0[2]:=c20:
end proc:

 

soln('parameters'=[1,0,0.1]);soln(0.1);

[alpha = 1., beta = 0., u = .1]

[t = .1, ca(t) = HFloat(0.9895549324188543), cb(t) = HFloat(0.009898024276129616)]

 

 

ssdsolve:=proc(x)
interface(warnlevel=0):

#if  type(x,Vector)#if type is not needed for this problem, might be needed for other problems
#then


local z1,n1,i,c10,c20,dt,u;
global soln,nvar;
dt:=evalf(1.0/nvar):
c10:=1.0:c20:=0.0:
for i from 1 to nvar do
u:=x[i]:
soln('parameters'=[c10,c20,u]):
z1:=soln(dt):
c10:=subs(z1,ca(t)):c20:=subs(z1,cb(t)):
od:
-c20;
 #else 'procname'(args):

#end if:

end proc:

 

 

ssRK2:=proc(x)#based on RK2
#interface(warnlevel=0):
option hfloat;
#if  type(x,Vector)
#then

local z1,n1,i,c10,c20,dt,u,NN,y0;
global nvar,RK2h;
y0:=Array(1..2,[1.0,0.0],datatype=float[8]):
#y0[1]:=1.0:y0[2]:=0.0:
dt:=evalf(1.0/nvar):
NN:=256*2/nvar;#NN is hardcode based on the assumption that nvar will be a multiple of 2 <=256


for i from 1 to nvar do
u:=x[i]:
evalhf(RK2h(NN,u,dt,y0)):
od:
-y0[2];
 #else 'procname'(args):

#end if:

end proc:

nvar:=2;

2

ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float):

bl := Vector(nvar,[seq(0.,i=1..nvar)],datatype=float):bu := Vector(nvar,[seq(5.,i=1..nvar)],datatype=float):

infolevel[Optimization]:=15;

15

CodeTools:-Usage(Optimization:-NLPSolve(nvar,evalhf(ssRK2),[],initialpoint=ic0,[bl,bu],optimalitytolerance=1e-6)):

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

NLPSolve: feasibility tolerance set to 0.1053671213e-7

NLPSolve: optimality tolerance set to 0.1e-5

NLPSolve: iteration limit set to 50

NLPSolve: infinite bound set to 0.10e21

NLPSolve: trying evalhf mode

NLPSolve: trying evalf mode

attemptsolution: number of major iterations taken 11

memory used=0.96MiB, alloc change=0 bytes, cpu time=16.00ms, real time=121.00ms, gc time=0ns

 

Calling the procedures based on RK2 with NLPSolve uses evalf for numerical gradient (fdiff). Providing a procedure for gradient solves this issue.

gradRK2 := proc (x,g)
option hfloat;
local base, i, xnew, del;
global nvar,ssRK2;
xnew:=Array(1..nvar,datatype=float[8]):
base := ssRK2(x);
for i to nvar do
xnew[i] := x[i]; end do;
for i to nvar do
del := max(1e-5,.1e-6*x[i]);
xnew[i] := xnew[i]+del;
g[i] := (ssRK2(xnew)-base)/del;
xnew[i] := xnew[i]-del;
end do;
end proc:

CodeTools:-Usage(Optimization:-NLPSolve(nvar,evalhf(ssRK2),[],initialpoint=ic0,[bl,bu],objectivegradient=gradRK2,optimalitytolerance=1e-6)):

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

NLPSolve: feasibility tolerance set to 0.1053671213e-7

NLPSolve: optimality tolerance set to 0.1e-5

NLPSolve: iteration limit set to 50

NLPSolve: infinite bound set to 0.10e21

NLPSolve: trying evalhf mode

attemptsolution: number of major iterations taken 11

memory used=184.77KiB, alloc change=0 bytes, cpu time=31.00ms, real time=110.00ms, gc time=0ns

Significant saving in memory used is seen. CPU time is also less, which is more apparent at larger valeus of nvar.

 

dsolvenumeric based codes work for optimization, but evalf is invoked possibly for both the objective and gradient

CodeTools:-Usage(Optimization:-NLPSolve(nvar,(ssdsolve),[],initialpoint=ic0,[bl,bu],optimalitytolerance=1e-6)):

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

NLPSolve: feasibility tolerance set to 0.1053671213e-7

NLPSolve: optimality tolerance set to 0.1e-5

NLPSolve: iteration limit set to 50

NLPSolve: infinite bound set to 0.10e21

NLPSolve: trying evalhf mode

NLPSolve: trying evalf mode

attemptsolution: number of major iterations taken 11

memory used=14.61MiB, alloc change=37.00MiB, cpu time=94.00ms, real time=213.00ms, gc time=46.88ms

 

Providing gradient procedure doesn't help with evalhf computaiton for dsolve/numeric based procedures.

graddsolve := proc (x,g)
local base, i, xnew, del;
global nvar,ssdsolve;
#xnew:=Vector(nvar,datatype=float[8]):
xnew:=Array(1..nvar,datatype=float[8]):
base := ssdsolve(x);
for i to nvar do
xnew[i] := x[i]; end do;
for i to nvar do
del := max(1e-5,.1e-6*x[i]);
xnew[i] := xnew[i]+del;
g[i] := (ssdsolve(xnew)-base)/del;
xnew[i] := xnew[i]-del;
end do;
end proc:

CodeTools:-Usage(Optimization:-NLPSolve(nvar,(ssdsolve),[],initialpoint=ic0,[bl,bu],objectivegradient=graddsolve,optimalitytolerance=1e-6)):

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

NLPSolve: feasibility tolerance set to 0.1053671213e-7

NLPSolve: optimality tolerance set to 0.1e-5

NLPSolve: iteration limit set to 50

NLPSolve: infinite bound set to 0.10e21

NLPSolve: trying evalhf mode

NLPSolve: trying evalf mode

attemptsolution: number of major iterations taken 11

memory used=3.82MiB, alloc change=0 bytes, cpu time=31.00ms, real time=129.00ms, gc time=0ns

 

Calling both RK2 and dsolve based procedures again to check the values and computation meterics

s1RK2:=CodeTools:-Usage(Optimization:-NLPSolve(nvar,evalhf(ssRK2),[],initialpoint=ic0,[bl,bu],objectivegradient=gradRK2,optimalitytolerance=1e-6)):

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

NLPSolve: feasibility tolerance set to 0.1053671213e-7

NLPSolve: optimality tolerance set to 0.1e-5

NLPSolve: iteration limit set to 50

NLPSolve: infinite bound set to 0.10e21

NLPSolve: trying evalhf mode

attemptsolution: number of major iterations taken 11

memory used=185.09KiB, alloc change=0 bytes, cpu time=16.00ms, real time=95.00ms, gc time=0ns

s1dsolve:=CodeTools:-Usage(Optimization:-NLPSolve(nvar,ssdsolve,[],initialpoint=ic0,[bl,bu],objectivegradient=graddsolve,optimalitytolerance=1e-6)):

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

NLPSolve: feasibility tolerance set to 0.1053671213e-7

NLPSolve: optimality tolerance set to 0.1e-5

NLPSolve: iteration limit set to 50

NLPSolve: infinite bound set to 0.10e21

NLPSolve: trying evalhf mode

NLPSolve: trying evalf mode

attemptsolution: number of major iterations taken 11

memory used=3.82MiB, alloc change=0 bytes, cpu time=31.00ms, real time=122.00ms, gc time=0ns

s1RK2[1];s1dsolve[1];

-.523900304316377463

-.523901163953022553

 

Next, a for loop is written to optimize for increasing values of nvar. One can see that evalhf is important for larger values of nvar. While dsolve/numeric is a superior code, not being able to use it in evalhf format is a significant weakness and should be addressed. Note that dsolve numeric evalutes procedures that are evaluated in evalhf or compiled form, so hopefully this is an easy fix.

infolevel[Optimization]:=0:

for j from 1 to 9 do
nvar:=2^(j-1):
ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float):
bl := Vector(nvar,[seq(0.,i=1..nvar)],datatype=float):bu := Vector(nvar,[seq(5.,i=1..nvar)],datatype=float):
soptRK[j]:=CodeTools:-Usage(Optimization:-NLPSolve(nvar,evalhf(ssRK2),[],initialpoint=ic0,[bl,bu],objectivegradient=gradRK2,optimalitytolerance=1e-6)):
print(2^(j-1),soptRK[j][1]);

od:

memory used=88.40KiB, alloc change=0 bytes, cpu time=0ns, real time=10.00ms, gc time=0ns

1, HFloat(-0.5008626988793192)

memory used=155.66KiB, alloc change=0 bytes, cpu time=16.00ms, real time=30.00ms, gc time=0ns

2, -.523900304316377463

memory used=175.36KiB, alloc change=0 bytes, cpu time=62.00ms, real time=80.00ms, gc time=0ns

4, -.535152497919956782

memory used=243.69KiB, alloc change=0 bytes, cpu time=156.00ms, real time=239.00ms, gc time=0ns

8, -.540546896131004706

memory used=260.09KiB, alloc change=0 bytes, cpu time=141.00ms, real time=260.00ms, gc time=0ns

16, -.542695734426874465

memory used=385.84KiB, alloc change=0 bytes, cpu time=313.00ms, real time=545.00ms, gc time=0ns

32, -.542932877726400531

memory used=0.65MiB, alloc change=0 bytes, cpu time=734.00ms, real time=1.18s, gc time=0ns

64, -.543011976841572652

memory used=1.24MiB, alloc change=0 bytes, cpu time=1.55s, real time=2.40s, gc time=0ns

128, -.543035276649319276

memory used=2.99MiB, alloc change=0 bytes, cpu time=3.45s, real time=5.92s, gc time=0ns

256, -.543041496228812814

for j from 1 to 6 do
nvar:=2^(j-1):
ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float):
bl := Vector(nvar,[seq(0.,i=1..nvar)],datatype=float):bu := Vector(nvar,[seq(5.,i=1..nvar)],datatype=float):
soptdsolve[j]:=CodeTools:-Usage(Optimization:-NLPSolve(nvar,evalf(ssdsolve),[],initialpoint=ic0,[bl,bu],objectivegradient=graddsolve,optimalitytolerance=1e-6)):
print(2^(j-1),soptdsolve[j][1]);

od:

memory used=0.66MiB, alloc change=0 bytes, cpu time=16.00ms, real time=11.00ms, gc time=0ns

1, HFloat(-0.5008631947224957)

memory used=3.79MiB, alloc change=0 bytes, cpu time=15.00ms, real time=52.00ms, gc time=0ns

2, -.523901163953022553

memory used=21.28MiB, alloc change=0 bytes, cpu time=78.00ms, real time=236.00ms, gc time=0ns

4, -.535153942647626391

memory used=144.82MiB, alloc change=-4.00MiB, cpu time=469.00ms, real time=1.60s, gc time=46.88ms

8, -.540549407239521607

memory used=347.30MiB, alloc change=16.00MiB, cpu time=1.27s, real time=3.45s, gc time=140.62ms

16, -.542699055038344258

memory used=1.33GiB, alloc change=-8.00MiB, cpu time=7.66s, real time=15.92s, gc time=750.00ms

32, -.542936165630524603

 

SS:=[seq(soptRK[j][1],j=1..9)];

[HFloat(-0.5008626988793192), -.523900304316377463, -.535152497919956782, -.540546896131004706, -.542695734426874465, -.542932877726400531, -.543011976841572652, -.543035276649319276, -.543041496228812814]

 

E1:=[seq(SS[i]-SS[i+1],i=1..nops(SS)-1)];

[HFloat(0.02303760543705824), 0.11252193603580e-1, 0.5394398211048e-2, 0.2148838295869e-2, 0.237143299527e-3, 0.79099115172e-4, 0.23299807746e-4, 0.6219579494e-5]

To get 6 Digits accuracy we need nvar = 256 which may not be attainable with dsolve/numeric approach unless we are able to call it evalhf format.

 

 


 

Download ImportanceofevalhfNLPSolve.mw

dsolve/numeric + NLPSolve shows inconsistent behavior. This combination is important for parameter estimation and optimal control. Can anyone fix this? Hopefully, I am not making a mistake.

restart:

 

Test code written by Dr. Venkat Subramanian at UT Austin, 05/31/2023. This code uses CVP approach (piecwise constant) to perform optimal control. NLPSolve combination with dsolve numeric parametric form is buggy and fails for some values of nvar, and works for some values of nvar. Ideally increasing nvar should show convergence with respect to the objective function.

restart:

Digits:=15;

15

 

eqodes:=[diff(ca(t),t)=-(u+u^2/2)*1.0*ca(t),diff(cb(t),t)=1.0*u*ca(t)-0.1*cb(t)];

[diff(ca(t), t) = -1.0*(u+(1/2)*u^2)*ca(t), diff(cb(t), t) = 1.0*u*ca(t)-.1*cb(t)]

soln:=dsolve({op(eqodes),ca(0)=alpha,cb(0)=beta},type=numeric,'parameters'=[alpha,beta,u],compile=true,savebinary=true):

 

 

ss:=proc(x)
interface(warnlevel=0):
#if  type(x[1],numeric)
if  type(x,Vector)
then

local z1,n1,i,c10,c20,dt,u;
global soln,nvar;
dt:=evalf(1.0/nvar):
c10:=1.0:c20:=0.0:
for i from 1 to nvar do
u:=x[i]:
soln('parameters'=[c10,c20,u]):
z1:=soln(dt):
c10:=subs(z1,ca(t)):c20:=subs(z1,cb(t)):
od:
-c20;
 else 'procname'(args):

end if:

end proc:

 

nvar:=3;#code works for nvar:=3, but not for nvar:=2

3

ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(3, {(1) = .1000000000000000, (2) = .1000000000000000, (3) = .1000000000000000})

 

ss(ic0);

HFloat(-0.09025793628011441)

bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(3, {(1) = 0., (2) = 0., (3) = 0.})

Vector[column](%id = 36893491136404053036)

infolevel[all]:=1;

1

Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 3

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

[-.531453381886523, Vector(3, {(1) = .8114345197312305, (2) = 1.189413022736622, (3) = 2.427689509710160})]

nvar:=2;#code works for nvar:=3, but not for nvar:=2

2

ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(2, {(1) = .1000000000000000, (2) = .1000000000000000})

 

ss(ic0);

HFloat(-0.09025793011810783)

bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(2, {(1) = 0., (2) = 0.})

Vector[column](%id = 36893491136437818540)

infolevel[all]:=1;

1

Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 2

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

Error, (in Optimization:-NLPSolve) no improved point could be found

 

nvar:=5;#code works for nvar:=3, but not for nvar:=2

5

ic0:=Vector(nvar,[seq(0.1,i=1..nvar)],datatype=float[8]);

Vector(5, {(1) = .1000000000000000, (2) = .1000000000000000, (3) = .1000000000000000, (4) = .1000000000000000, (5) = .1000000000000000})

 

ss(ic0);

HFloat(-0.09025792639212991)

bl := Vector(nvar,[seq(0.,i=1..nvar)]);bu := Vector(nvar,[seq(5.,i=1..nvar)]);

Vector(5, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0.})

Vector[column](%id = 36893491136472713804)

infolevel[all]:=1;

1

Optimization:-NLPSolve(nvar,ss,[],initialpoint=ic0,[bl,bu]);

NLPSolve: calling NLP solver

NLPSolve: using method=sqp

NLPSolve: number of problem variables 5

NLPSolve: number of nonlinear inequality constraints 0

NLPSolve: number of nonlinear equality constraints 0

NLPSolve: number of general linear constraints 0

[-.537338871783244, Vector(5, {(1) = .7435767224059456, (2) = .9013440906589676, (3) = 1.148772394841535, (4) = 1.629079877401040, (5) = 3.222801229872320})]

 


 

Download test.mw

A robust largescale index-1 DAE (and stiff ODE) solver has been developed in Maple.

PS: These codes run only with the Java version of Maple. Codes are saved as *.mws for the ease of typing. One can download the files and save them as *.mw as well.
Details about our approach are provided in a paper just submitted. 
See arxiv at
https://arxiv.org/abs/2212.02630

We encourage everyone to test these codes and report bugs. All the examples can be run with Intel's Pardiso (provided users have access to libraries) by calling DAESolverP.txt and by calling "IMPDAEP" instead of "IMPDAE". This is useful for large-scale problems. The symbolic capability and ListTools search capability of Maple are very good and can be used for developing optimization solvers as well.

I would like to know if CPU time/memory usage can be reduced significantly. In particular,  for examples 5 and 6. Some ways to contribute include

(1) Running the code in evalhf or compiled form. This may be hard.
(2) Providing options to run other parallel open-source linear solvers (eg., MUMPS).

(3) Other examples that show the use of the developed solver. We are able to solve > 100,000 DAEs.

(4) Helping in converting the code to Maple 14 or earlier (by doing sparse LU Decomposition. Just using LinearSolve will slow down the code).

 

Please avoid ~,*, etc (shortcuts) unless it improves the speed of calculation.
Thanks 
Dr. Venkat Subramanian

Download examples_1_3.mws

Download example_4.mws

Download example_5.mws

Download example_6.mws

 

mw format files are given below
Download examples_1_3.mw

Download example_4.mw

Download example_5.mw

Download example_6.mw

DAESolver.txt

DAESolverP.txt

Details about our approach for phase-field models are provided in a paper just submitted. 
https://ecsarxiv.org/k2vu6/
 

Maple code (based on UMFPACK linearsolver and different compiled procedures for marching in time) is given here. Run the code for small values of NN and MM and keep increasing them. The goal is to do 100, 200, 400, etc and still get the code to run very fast.
 

restart:

t11:=time():

t12:=time[real]():

Digits:=15;

15

(1)

NN is the number of node points (elements) in the X and MM is the number of elements in the Y direction. delta is the applied current density at the top (Y =1). tf is the final time for simulation. vel is the velocity constant v in the paper. ki0 is the scaled exchange current density k in the paper. This code can be run for positive values of delta. This simulates plating. At the end of simulation, changing delta to negative values and rerunning the code will automatically used the geometry at the end of plating.
Ydatstore stores the geometry at every point in time. Phiaveadd stores the total liquid phase in the domain at any point in time.
Users can change NN, delta, tf, vel, ki0, MM just in this line and choose Edit execute worksheet to run for different design parameters.

Users can modify the call for y0proc for choosing different models.

Users can modify the call for HF to run first-order upwind, ENO2 or WENO3 methods. NN and MM should be even numbers.

NN:=100;delta:=0.1;tf:=2.0;vel:=1.0;ki0:=1.0;MM:=NN;

100

 

.1

 

2.0

 

1.0

 

1.0

 

100

(2)

N:=NN+2:h:=1.0/NN:M:=MM+2;Ny:=MM;

102

 

100

(3)

 

Initial geometry, Model 1, semicircle in a square

y0proc1:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;

N:=NN+2:h:=1.0/NN:w:=h/2:

M:=MM+2;Ny:=MM;

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=xx^2+yy^2;
Y00[i,j]:=max(1e-9,0.5+0.5*tanh((sqrt(rr)-0.3)/w/sqrt(2.0))):
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

 

Model 2, square in a square

y0proc2:=proc(NN,MM,Y00)# square inside a circle, Model 2
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;

N:=NN+2:h:=1.0/NN:w:=h/2:

M:=MM+2;Ny:=MM;

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
if xx <=0.3 and yy<=0.3 then Y00[i,j]:=1e-9; else Y00[i,j]:=1.0:end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

 

 

Initial geometry, Model 3, electrodeposition problem trenches and via

y0proc3:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;

N:=NN+2:h:=1.0/NN:w:=h/2:

M:=MM+2;Ny:=MM;

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
if abs(xx-0.5) >0.2 and yy<=0.5 then Y00[i,j]:=1e-9; else Y00[i,j]:=1.0:end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

 

 

Initial geometry, Model 4, Gaussian Seed at the bottom

y0proc4:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;

N:=NN+2:h:=1.0/NN:w:=h/2:

M:=MM+2;Ny:=MM;

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=0.1+0.1*exp(-500.*(xx-0.5)^2);
Y00[i,j]:=max(1e-9,0.5+0.5*tanh((yy-rr)/w/sqrt(2.0))):
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

 

y0proc:=eval(y0proc1):#choose different models using y0proc2, etc.

Y00:=Matrix(1..N,1..M,datatype=float[8]):

evalhf(y0proc(NN,MM,Y00)):

if delta<0 then read("Y0data.m"):end:

if delta<0 then read("tdata.m"):end:

tf;

2.0

(4)

p0:=plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]):

printlevel:=2;

2

(5)

eq:=Array(1..N,1..M):

Nx:=NN:

 

 

Next, boundary conditiosn at X  = 0, X =1, Y = 0, Y = 1 are specified below, but these equations are not used and optimally coded inside the procedure Eqs11.

for i from 1 to 1 do for j from 1 to M do eq[i,j]:=-Phi[i,j]+Phi[i+1,j]:od:od:

for i from N to N do for j from 1 to M do eq[i,j]:=Phi[i-1,j]-Phi[i,j]:od:od:

for i from 1 to N do for j from 1 to 1 do eq[i,j]:=Phi[i,j+1]-Phi[i,j]:od:od:

for i from 1 to N do for j from M to M do eq[i,j]:=Phi[i,j-1]-Phi[i,j]+delta*h:od:od:

Eqs:=Vector(N*(M)):

N,M;

102, 102

(6)

 

Residues at different points in X and Y are coded in the Eqs11 procedure. Y0 is the input phase-field parameter (2D Matrix). The potential y is a vector to expedite the calculation of residue.

Eqs11:=proc(Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,ff::Vector(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
global N,M;
h:=1.0/(N-2):

for i from 2 to N-1 do for j from 2 to M-1 do
i1:=i+(j-1)*N:
ff[i1]:=
(Y0[i,j]+Y0[i,j+1])*(y[i1+N]-y[i1])
-(Y0[i,j]+Y0[i,j-1])*(y[i1]-y[i1-N])
+(Y0[i+1,j]+Y0[i,j])*(y[i1+1]-y[i1])
-(Y0[i,j]+Y0[i-1,j])*(y[i1]-y[i1-1])
-ki0*y[i1]*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
od:od:
for i from 1 to 1 do for j from 1 to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+1]:
od:od:
for i from N to N do for j from 1 to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-1]:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+N]:
od:od:
for i from 1 to N do for j from M to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-N]+delta*h:
od:od:
#print(1);
#ff;
end proc:

 

The Jacobian for the residues is coded int he procedure Jac1.

Jac1:=proc(Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,j00::Matrix(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
global N,M;
h:=1.0/(N-2):

for i from 2 to N-1 do for j from 2 to M-1 do
i1:=i+(j-1)*N:
j00[i1,i1]:=-4*Y0[i,j]-Y0[i,j+1]-Y0[i,j-1]-Y0[i+1,j]-Y0[i-1,j]-ki0*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
j00[i1,i1+1]:=Y0[i+1,j]+Y0[i,j]:j00[i1,i1-1]:=Y0[i-1,j]+Y0[i,j]:
j00[i1,i1+N]:=Y0[i,j+1]+Y0[i,j]:j00[i1,i1-N]:=Y0[i,j-1]+Y0[i,j]:
od:od:
for i from 1 to 1 do for j from 1 to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1+1]:=1.00:
#ff[i1]:=-y[i1]+y[i1+1]:
od:od:
for i from N to N do for j from 1 to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1-1]:=1.00:
#ff[i1]:=-y[i1]+y[i1-1]:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*N:
#ff[i1]:=-y[i1]+y[i1+N]:
j00[i1,i1]:=-1.0:j00[i1,i1+N]:=1.00:
od:od:
for i from 1 to N do for j from M to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1-N]:=1.00:
#ff[i1]:=-y[i1]+y[i1-N]+delta*h:
od:od:
#print(1);
#ff;
end proc:

Eqs11:=Compiler:-Compile(Eqs11):

 

ntot:=N*(M);

10404

(7)

Ntot:=ntot;

10404

(8)

 

First-order Upwind method is coded in the procedure UpW.

UpW:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8];
#e1:=1e-15:

h:=1.0/(N-2):

vv0:=v0:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:


vx1:=(Y00[i,j]-Y00[i-1,j])/h:
vx2:=(Y00[i+1,j]-Y00[i,j])/h:


vy1:=(Y00[i,j]-Y00[i,j-1])/h:
vy2:=(Y00[i,j+1]-Y00[i,j])/h:

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
F0[i,j]:=nx:
od:od:

end proc:

 

Second-order ENO2 method is coded in the procedure ENO2.

ENO2:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],alpha::float[8],beta::float[8];


h:=1.0/(N-2):
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:

sdx:=(Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j])/h:sdy:=(Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1])/h:
vxb:=0:vxf:=0:vyb:=0:vyf:=0:
if i = 2 then sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j])/h:else sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j])/h:end:
if sdx*sdxb>=0 then s1x:=1.0 else s1x:=0.0:end:
vx1:=(Y00[i,j]-Y00[i-1,j])/h+0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxb)):

if i = N-1 then sdxf:=(Y00[i+1,j]-2*Y00[i+1,j]+Y00[i,j])/h:else sdxf:=(Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j])/h:end:
if sdx*sdxf>=0 then s1x:=1.0 else s1x:=0.0:end:
vx2:=(Y00[i+1,j]-Y00[i,j])/h-0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxf)):


if j = 2 then sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1])/h:else sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2])/h:end:
if sdy*sdyb>=0 then s1y:=1.0 else s1y:=0.0:end:
vy1:=(Y00[i,j]-Y00[i,j-1])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyb)):

if j = M-1 then sdyf:=(Y00[i,j+1]-2*Y00[i,j+1]+Y00[i,j])/h:else sdyf:=(Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j])/h:end:
if sdy*sdyf>=0 then s1y:=1.0 else s1y:=0.0:end:
vy2:=(Y00[i,j+1]-Y00[i,j])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyf)):

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):

F0[i,j]:=nx:
od:od:

end proc:

 

Third-order WENO3 method is coded in the procedure WENO3.

WENO3:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8],e1::float[8];
e1:=1e-6:

h:=1.0/(N-2):

vv0:=v0:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:
phix:=(Y00[i+1,j]-Y00[i-1,j])/2/h:
phiy:=(Y00[i,j+1]-Y00[i,j-1])/2/h:

if i = 2 then sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j]: else sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j]:end:
sd:=Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j]:
if i = N-1 then sdf:=Y00[i,j]-2*Y00[i+1,j]+Y00[i+1,j]: else sdf:=Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):
vx1:=phix-0.5*w1/h*(sd-sdb):
vx2:=phix-0.5*w2/h*(sdf-sd):

if j = 2 then sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1]:else sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2]:end:
sd:=Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1]:
if j = M-1 then sdf:=Y00[i,j]-2*Y00[i,j+1]+Y00[i,j+1]: else sdf:=Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):

vy1:=phiy-0.5*w1/h*(sd-sdb):
vy2:=phiy-0.5*w2/h*(sdf-sd):

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
#nx:=sqrt(vx1^2+vx2^2+vy1^2+vy^2):

F0[i,j]:=nx:
od:od:

end proc:

 

UpW:=Compiler:-Compile(UpW):

ENO2:=Compiler:-Compile(ENO2):

WENO3:=Compiler:-Compile(WENO3):

F0:=Matrix(1..N,1..M,datatype=float[8]):

h/0.1/2;

0.500000000000000e-1

(9)

 

Nx;

100

(10)

PhiAdd:=proc(N::integer,Phi0::Vector(datatype=float[8]),db::Vector(datatype=float[8]))
local i::integer;
for i from 1 to N do Phi0[i]:=Phi0[i]+db[i]:od:
end proc:

PhiAdd:=Compiler:-Compile(PhiAdd):

Ntot:=ntot:

Phi0:=Vector(1..Ntot,datatype=float[8]):

printlevel:=1:

ff:=copy(Phi0):

 

Phi0:=Vector(1..Ntot,datatype=float[8]):

j00:=Matrix(1..Ntot,1..Ntot,datatype=float[8],storage=sparse):

evalf(Eqs11(Y00,Phi0,delta,ki0,ff));

0.100000000000000e-2

(11)

Jac1(Y00,Phi0,delta,ki0,j00):

db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):

Phi0:=Phi0+db:

V[0]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*delta;

HFloat(0.31581237885872887)

(12)

TT[0]:=0;tt:=0:

0

(13)

vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):

dt:=min(h/vv0/vel/ki0,tf-tt);

HFloat(0.031202770994846224)

(14)

Ymid:=Matrix(1..N,1..M,datatype=float[8]):

phiaveAdd:=proc(Y00::Matrix(datatype=float[8]),N::integer,M::integer,Ny::integer)
local i::integer,j::integer,phiave::float;
phiave:=0.0:
for i from 2 to N-1 do for j from 2 to M-1 do phiave:=phiave+Y00[i,j]:od:od:
phiave/(N-2)/(M-2);
end proc:

phiaveAdd:=Compiler:-Compile(phiaveAdd):

Ny;M;

100

 

102

(15)

phiaveAdd(Y00,N,M,Ny);

.929280193991299241

(16)

 

The three stages of  SSR-RK3 are stored in EFAdd, EFAdd2 and EFAdd3

EFAdd:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
#for i from 2 to N-1 do for j from 2 to M-1 do Ymid[i,j]:=Y00[i,j]-dt*F0[i,j]*Phi0[i+(j-1)*N]:od:od:
for i from 2 to N-1 do for j from 2 to M-1 do Ymid[i,j]:=max(1e-9,Y00[i,j]-dt*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Ymid[i,1]:=Ymid[i,2]:
Ymid[i,M]:=Ymid[i,M-1]:
od:
for j from 1 to M do
Ymid[1,j]:=Ymid[2,j]:
Ymid[N,j]:=Ymid[N-1,j]:
od:
end proc:

EFAdd:=Compiler:-Compile(EFAdd):

EFAdd2:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Ymid[i,j]:=max(1e-9,Y00[i,j]*3/4.+Ymid[i,j]/4.-dt/4.*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Ymid[i,1]:=Ymid[i,2]:
Ymid[i,M]:=Ymid[i,M-1]:
od:
for j from 1 to M do
Ymid[1,j]:=Ymid[2,j]:
Ymid[N,j]:=Ymid[N-1,j]:
od:
end proc:

EFAdd2:=Compiler:-Compile(EFAdd2):

EFAdd3:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Y00[i,j]:=max(1e-9,Y00[i,j]*1/3.+Ymid[i,j]*2/3.-dt*2/3.*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

EFAdd3:=Compiler:-Compile(EFAdd3):

 

YdatStore:=proc(Y00::Matrix(datatype=float[8]),Ydat::Array(datatype=float[8]),N::integer,M::integer,jj::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Ydat[jj,i,j]:=Y00[i,j]:od:od:
end proc:

YdatStore:=Compiler:-Compile(YdatStore):

Phiave[0]:=phiaveAdd(Y00,N,M,Ny);

.929280193991299241

(17)

 

Nt:=round(tf/dt)+50;

114

(18)

printlevel:=1:

Ymid:=copy(Y00):

Ydat:=Array(1..Nt+1,1..N,1..M,datatype=float[8]):

YdatStore(Y00,Ydat,N,M,1);

1.

(19)

 

ii:=0:TT[0];tt;

0

 

0

(20)

 

Different upwind schemes can be called by assign WENO3, UpW or ENO3 scheme.

HF:=eval(WENO3):

#HF:=eval(UpW):

#HF:=eval(ENO2):

 

A while loop is written from t=0 to t= tf.

while tt<tf do
HF(Y00,Phi0,F0,evalf(dt),N,M,delta):
EFAdd(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Ymid,Phi0,delta,ki0,ff));
Jac1(Ymid,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd2(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Ymid,Phi0,delta,ki0,ff));
Jac1(Ymid,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd3(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Y00,Phi0,delta,ki0,ff));
Jac1(Y00,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
ii:=ii+1:#print(i);
V[ii]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*delta;#print(ii,V[ii]);
TT[ii]:=TT[ii-1]+dt:tt:=tt+dt:
YdatStore(Y00,Ydat,N,M,ii+1);
Phiave[ii]:=phiaveAdd(Y00,N,M,Ny);
vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
dt:=min(h/vv0/vel/ki0,tf-tt);#print(tt,ii,dt);
#YdatStore(Y00,Ydat,N,M,i+1);
#Phiave[i]:=phiaveAdd(Y00,N,M,Ny);
end:

 

nt:=ii;

46

(21)

V[nt];

HFloat(0.17021428328406052)

(22)

 

Voltage time curves are plotted below. Voltage is measured at X = 0.5, Y = 1.

plot([[seq([TT[ii],V[ii]],ii=0..nt)]],style=point);

 

 

Voltage at the end of plating, cpu time can be documented as

V[nt];time[real]()-t12;time()-t11;NN;

HFloat(0.17021428328406052)

 

19.042

 

14.719

 

100

(23)

 

p1:=plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]):

tf:=TT[nt];

HFloat(2.0)

(24)

Contour plots at t= 0 and t = 2.0 (at the of plating are given below)

plots:-display({p0});plots:-display({p1});

 

 

 

 

The liquid phase content as a funciton of time is plotted below

plot([[seq([TT[ii],Phiave[ii]],ii=0..nt)]],style=point);

 

 

save Y00,"Y0data.m";

save tf,"tdata.m";

read("Y0data.m"):

plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]);

 

 


 

Download Phasefield2DBaseCodeParametric.mw
 

restart:

t11:=time():

t12:=time[real]():

Digits:=15;

15

(1)

NN is the number of node points (elements) in the X and MM is the number of elements in the Y direction. delta is the applied current density at the top (Y =1). tf is the final time for simulation. vel is the velocity constant v in the paper. ki0 is the scaled exchange current density k in the paper. This code can be run for positive values of delta. This simulates plating. At the end of simulation, changing delta to negative values and rerunning the code will automatically used the geometry at the end of plating.
Ydatstore stores the geometry at every point in time. Phiaveadd stores the total liquid phase in the domain at any point in time.
Users can change NN, delta, tf, vel, ki0, MM just in this line and choose Edit execute worksheet to run for different design parameters.

Users can modify the call for y0proc for choosing different models.

Users can modify the call for HF to run first-order upwind, ENO2 or WENO3 methods. NN and MM should be even numbers.

NN:=100;delta:=0.1;tf:=2.0;vel:=1.0;ki0:=1.0;MM:=NN;

100

 

.1

 

2.0

 

1.0

 

1.0

 

100

(2)

N:=NN+2:h:=1.0/NN:M:=MM+2;Ny:=MM;

102

 

100

(3)

 

Initial geometry, Model 1, semicircle in a square

y0proc1:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;

N:=NN+2:h:=1.0/NN:w:=h/2:

M:=MM+2;Ny:=MM;

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=xx^2+yy^2;
Y00[i,j]:=max(1e-9,0.5+0.5*tanh((sqrt(rr)-0.3)/w/sqrt(2.0))):
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

 

Model 2, square in a square

y0proc2:=proc(NN,MM,Y00)# square inside a circle, Model 2
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;

N:=NN+2:h:=1.0/NN:w:=h/2:

M:=MM+2;Ny:=MM;

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
if xx <=0.3 and yy<=0.3 then Y00[i,j]:=1e-9; else Y00[i,j]:=1.0:end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

 

 

Initial geometry, Model 3, electrodeposition problem trenches and via

y0proc3:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;

N:=NN+2:h:=1.0/NN:w:=h/2:

M:=MM+2;Ny:=MM;

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
if abs(xx-0.5) >0.2 and yy<=0.5 then Y00[i,j]:=1e-9; else Y00[i,j]:=1.0:end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

 

 

Initial geometry, Model 4, Gaussian Seed at the bottom

y0proc4:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;

N:=NN+2:h:=1.0/NN:w:=h/2:

M:=MM+2;Ny:=MM;

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=0.1+0.1*exp(-500.*(xx-0.5)^2);
Y00[i,j]:=max(1e-9,0.5+0.5*tanh((yy-rr)/w/sqrt(2.0))):
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

 

y0proc:=eval(y0proc1):#choose different models using y0proc2, etc.

Y00:=Matrix(1..N,1..M,datatype=float[8]):

evalhf(y0proc(NN,MM,Y00)):

if delta<0 then read("Y0data.m"):end:

if delta<0 then read("tdata.m"):end:

tf;

2.0

(4)

p0:=plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]):

printlevel:=2;

2

(5)

eq:=Array(1..N,1..M):

Nx:=NN:

 

 

Next, boundary conditiosn at X  = 0, X =1, Y = 0, Y = 1 are specified below, but these equations are not used and optimally coded inside the procedure Eqs11.

for i from 1 to 1 do for j from 1 to M do eq[i,j]:=-Phi[i,j]+Phi[i+1,j]:od:od:

for i from N to N do for j from 1 to M do eq[i,j]:=Phi[i-1,j]-Phi[i,j]:od:od:

for i from 1 to N do for j from 1 to 1 do eq[i,j]:=Phi[i,j+1]-Phi[i,j]:od:od:

for i from 1 to N do for j from M to M do eq[i,j]:=Phi[i,j-1]-Phi[i,j]+delta*h:od:od:

Eqs:=Vector(N*(M)):

N,M;

102, 102

(6)

 

Residues at different points in X and Y are coded in the Eqs11 procedure. Y0 is the input phase-field parameter (2D Matrix). The potential y is a vector to expedite the calculation of residue.

Eqs11:=proc(Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,ff::Vector(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
global N,M;
h:=1.0/(N-2):

for i from 2 to N-1 do for j from 2 to M-1 do
i1:=i+(j-1)*N:
ff[i1]:=
(Y0[i,j]+Y0[i,j+1])*(y[i1+N]-y[i1])
-(Y0[i,j]+Y0[i,j-1])*(y[i1]-y[i1-N])
+(Y0[i+1,j]+Y0[i,j])*(y[i1+1]-y[i1])
-(Y0[i,j]+Y0[i-1,j])*(y[i1]-y[i1-1])
-ki0*y[i1]*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
od:od:
for i from 1 to 1 do for j from 1 to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+1]:
od:od:
for i from N to N do for j from 1 to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-1]:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+N]:
od:od:
for i from 1 to N do for j from M to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-N]+delta*h:
od:od:
#print(1);
#ff;
end proc:

 

The Jacobian for the residues is coded int he procedure Jac1.

Jac1:=proc(Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,j00::Matrix(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
global N,M;
h:=1.0/(N-2):

for i from 2 to N-1 do for j from 2 to M-1 do
i1:=i+(j-1)*N:
j00[i1,i1]:=-4*Y0[i,j]-Y0[i,j+1]-Y0[i,j-1]-Y0[i+1,j]-Y0[i-1,j]-ki0*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
j00[i1,i1+1]:=Y0[i+1,j]+Y0[i,j]:j00[i1,i1-1]:=Y0[i-1,j]+Y0[i,j]:
j00[i1,i1+N]:=Y0[i,j+1]+Y0[i,j]:j00[i1,i1-N]:=Y0[i,j-1]+Y0[i,j]:
od:od:
for i from 1 to 1 do for j from 1 to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1+1]:=1.00:
#ff[i1]:=-y[i1]+y[i1+1]:
od:od:
for i from N to N do for j from 1 to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1-1]:=1.00:
#ff[i1]:=-y[i1]+y[i1-1]:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*N:
#ff[i1]:=-y[i1]+y[i1+N]:
j00[i1,i1]:=-1.0:j00[i1,i1+N]:=1.00:
od:od:
for i from 1 to N do for j from M to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1-N]:=1.00:
#ff[i1]:=-y[i1]+y[i1-N]+delta*h:
od:od:
#print(1);
#ff;
end proc:

Eqs11:=Compiler:-Compile(Eqs11):

 

ntot:=N*(M);

10404

(7)

Ntot:=ntot;

10404

(8)

 

First-order Upwind method is coded in the procedure UpW.

UpW:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8];
#e1:=1e-15:

h:=1.0/(N-2):

vv0:=v0:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:


vx1:=(Y00[i,j]-Y00[i-1,j])/h:
vx2:=(Y00[i+1,j]-Y00[i,j])/h:


vy1:=(Y00[i,j]-Y00[i,j-1])/h:
vy2:=(Y00[i,j+1]-Y00[i,j])/h:

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
F0[i,j]:=nx:
od:od:

end proc:

 

Second-order ENO2 method is coded in the procedure ENO2.

ENO2:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],alpha::float[8],beta::float[8];


h:=1.0/(N-2):
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:

sdx:=(Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j])/h:sdy:=(Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1])/h:
vxb:=0:vxf:=0:vyb:=0:vyf:=0:
if i = 2 then sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j])/h:else sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j])/h:end:
if sdx*sdxb>=0 then s1x:=1.0 else s1x:=0.0:end:
vx1:=(Y00[i,j]-Y00[i-1,j])/h+0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxb)):

if i = N-1 then sdxf:=(Y00[i+1,j]-2*Y00[i+1,j]+Y00[i,j])/h:else sdxf:=(Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j])/h:end:
if sdx*sdxf>=0 then s1x:=1.0 else s1x:=0.0:end:
vx2:=(Y00[i+1,j]-Y00[i,j])/h-0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxf)):


if j = 2 then sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1])/h:else sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2])/h:end:
if sdy*sdyb>=0 then s1y:=1.0 else s1y:=0.0:end:
vy1:=(Y00[i,j]-Y00[i,j-1])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyb)):

if j = M-1 then sdyf:=(Y00[i,j+1]-2*Y00[i,j+1]+Y00[i,j])/h:else sdyf:=(Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j])/h:end:
if sdy*sdyf>=0 then s1y:=1.0 else s1y:=0.0:end:
vy2:=(Y00[i,j+1]-Y00[i,j])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyf)):

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):

F0[i,j]:=nx:
od:od:

end proc:

 

Third-order WENO3 method is coded in the procedure WENO3.

WENO3:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8],e1::float[8];
e1:=1e-6:

h:=1.0/(N-2):

vv0:=v0:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:
phix:=(Y00[i+1,j]-Y00[i-1,j])/2/h:
phiy:=(Y00[i,j+1]-Y00[i,j-1])/2/h:

if i = 2 then sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j]: else sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j]:end:
sd:=Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j]:
if i = N-1 then sdf:=Y00[i,j]-2*Y00[i+1,j]+Y00[i+1,j]: else sdf:=Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):
vx1:=phix-0.5*w1/h*(sd-sdb):
vx2:=phix-0.5*w2/h*(sdf-sd):

if j = 2 then sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1]:else sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2]:end:
sd:=Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1]:
if j = M-1 then sdf:=Y00[i,j]-2*Y00[i,j+1]+Y00[i,j+1]: else sdf:=Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):

vy1:=phiy-0.5*w1/h*(sd-sdb):
vy2:=phiy-0.5*w2/h*(sdf-sd):

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
#nx:=sqrt(vx1^2+vx2^2+vy1^2+vy^2):

F0[i,j]:=nx:
od:od:

end proc:

 

UpW:=Compiler:-Compile(UpW):

ENO2:=Compiler:-Compile(ENO2):

WENO3:=Compiler:-Compile(WENO3):

F0:=Matrix(1..N,1..M,datatype=float[8]):

h/0.1/2;

0.500000000000000e-1

(9)

 

Nx;

100

(10)

PhiAdd:=proc(N::integer,Phi0::Vector(datatype=float[8]),db::Vector(datatype=float[8]))
local i::integer;
for i from 1 to N do Phi0[i]:=Phi0[i]+db[i]:od:
end proc:

PhiAdd:=Compiler:-Compile(PhiAdd):

Ntot:=ntot:

Phi0:=Vector(1..Ntot,datatype=float[8]):

printlevel:=1:

ff:=copy(Phi0):

 

Phi0:=Vector(1..Ntot,datatype=float[8]):

j00:=Matrix(1..Ntot,1..Ntot,datatype=float[8],storage=sparse):

evalf(Eqs11(Y00,Phi0,delta,ki0,ff));

0.100000000000000e-2

(11)

Jac1(Y00,Phi0,delta,ki0,j00):

db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):

Phi0:=Phi0+db:

V[0]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*delta;

HFloat(0.31581237885872887)

(12)

TT[0]:=0;tt:=0:

0

(13)

vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):

dt:=min(h/vv0/vel/ki0,tf-tt);

HFloat(0.031202770994846224)

(14)

Ymid:=Matrix(1..N,1..M,datatype=float[8]):

phiaveAdd:=proc(Y00::Matrix(datatype=float[8]),N::integer,M::integer,Ny::integer)
local i::integer,j::integer,phiave::float;
phiave:=0.0:
for i from 2 to N-1 do for j from 2 to M-1 do phiave:=phiave+Y00[i,j]:od:od:
phiave/(N-2)/(M-2);
end proc:

phiaveAdd:=Compiler:-Compile(phiaveAdd):

Ny;M;

100

 

102

(15)

phiaveAdd(Y00,N,M,Ny);

.929280193991299241

(16)

 

The three stages of  SSR-RK3 are stored in EFAdd, EFAdd2 and EFAdd3

EFAdd:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
#for i from 2 to N-1 do for j from 2 to M-1 do Ymid[i,j]:=Y00[i,j]-dt*F0[i,j]*Phi0[i+(j-1)*N]:od:od:
for i from 2 to N-1 do for j from 2 to M-1 do Ymid[i,j]:=max(1e-9,Y00[i,j]-dt*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Ymid[i,1]:=Ymid[i,2]:
Ymid[i,M]:=Ymid[i,M-1]:
od:
for j from 1 to M do
Ymid[1,j]:=Ymid[2,j]:
Ymid[N,j]:=Ymid[N-1,j]:
od:
end proc:

EFAdd:=Compiler:-Compile(EFAdd):

EFAdd2:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Ymid[i,j]:=max(1e-9,Y00[i,j]*3/4.+Ymid[i,j]/4.-dt/4.*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Ymid[i,1]:=Ymid[i,2]:
Ymid[i,M]:=Ymid[i,M-1]:
od:
for j from 1 to M do
Ymid[1,j]:=Ymid[2,j]:
Ymid[N,j]:=Ymid[N-1,j]:
od:
end proc:

EFAdd2:=Compiler:-Compile(EFAdd2):

EFAdd3:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Y00[i,j]:=max(1e-9,Y00[i,j]*1/3.+Ymid[i,j]*2/3.-dt*2/3.*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:

EFAdd3:=Compiler:-Compile(EFAdd3):

 

YdatStore:=proc(Y00::Matrix(datatype=float[8]),Ydat::Array(datatype=float[8]),N::integer,M::integer,jj::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Ydat[jj,i,j]:=Y00[i,j]:od:od:
end proc:

YdatStore:=Compiler:-Compile(YdatStore):

Phiave[0]:=phiaveAdd(Y00,N,M,Ny);

.929280193991299241

(17)

 

Nt:=round(tf/dt)+50;

114

(18)

printlevel:=1:

Ymid:=copy(Y00):

Ydat:=Array(1..Nt+1,1..N,1..M,datatype=float[8]):

YdatStore(Y00,Ydat,N,M,1);

1.

(19)

 

ii:=0:TT[0];tt;

0

 

0

(20)

 

Different upwind schemes can be called by assign WENO3, UpW or ENO3 scheme.

HF:=eval(WENO3):

#HF:=eval(UpW):

#HF:=eval(ENO2):

 

A while loop is written from t=0 to t= tf.

while tt<tf do
HF(Y00,Phi0,F0,evalf(dt),N,M,delta):
EFAdd(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Ymid,Phi0,delta,ki0,ff));
Jac1(Ymid,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd2(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Ymid,Phi0,delta,ki0,ff));
Jac1(Ymid,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd3(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Y00,Phi0,delta,ki0,ff));
Jac1(Y00,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
ii:=ii+1:#print(i);
V[ii]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*delta;#print(ii,V[ii]);
TT[ii]:=TT[ii-1]+dt:tt:=tt+dt:
YdatStore(Y00,Ydat,N,M,ii+1);
Phiave[ii]:=phiaveAdd(Y00,N,M,Ny);
vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
dt:=min(h/vv0/vel/ki0,tf-tt);#print(tt,ii,dt);
#YdatStore(Y00,Ydat,N,M,i+1);
#Phiave[i]:=phiaveAdd(Y00,N,M,Ny);
end:

 

nt:=ii;

46

(21)

V[nt];

HFloat(0.17021428328406052)

(22)

 

Voltage time curves are plotted below. Voltage is measured at X = 0.5, Y = 1.

plot([[seq([TT[ii],V[ii]],ii=0..nt)]],style=point);

 

 

Voltage at the end of plating, cpu time can be documented as

V[nt];time[real]()-t12;time()-t11;NN;

HFloat(0.17021428328406052)

 

19.042

 

14.719

 

100

(23)

 

p1:=plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]):

tf:=TT[nt];

HFloat(2.0)

(24)

Contour plots at t= 0 and t = 2.0 (at the of plating are given below)

plots:-display({p0});plots:-display({p1});

 

 

 

 

The liquid phase content as a funciton of time is plotted below

plot([[seq([TT[ii],Phiave[ii]],ii=0..nt)]],style=point);

 

 

save Y00,"Y0data.m";

save tf,"tdata.m";

read("Y0data.m"):

plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]);

 

 


 

Download Phasefield2DBaseCodeParametric.mw

 

1 2 3 Page 1 of 3