- Questions - MaplePrimes

Inconsistent behavior with dsolve/numeri...

Asked: Dr. Venkat Subramanian 351

May 31 2023

3 23

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;

>

>	eqodes:=[diff(ca(t),t)=-(u+u^2/2)1.0ca(t),diff(cb(t),t)=1.0uca(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

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

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

>

ss(ic0);

>	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;

>	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

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

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

>

ss(ic0);

>	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;

>	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

>	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);

>	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;

>	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

Will anyone be able to speed up this cod...

Asked: Dr. Venkat Subramanian 351 Product: Maple

December 07 2022

0 3

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

Will anyone be able to speed up this cod...

Asked: Dr. Venkat Subramanian 351 Product: Maple

April 22 2022

1 15

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;

(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;

(2)

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

(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;

(4)

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

>	printlevel:=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;

(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);

(7)

>	Ntot:=ntot;

(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;

(9)

>

Nx;

(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));

(11)

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

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

>	Phi0:=Phi0+db:

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

(12)

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

(13)

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

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

(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;

(15)

>	phiaveAdd(Y00,N,M,Ny);

(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);

(17)

>

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

(18)

>	printlevel:=1:

>	Ymid:=copy(Y00):

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

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

(19)

>

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

(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;

(21)

>

V[nt];

(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;

(23)

>

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

>	tf:=TT[nt];

(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..NNh+h/2,-h/2..MMh+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]);

>

Download Phasefield2DBaseCodeParametric.mw

>

restart:

>	t11:=time():

>	t12:=time[real]():

>	Digits:=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;

(2)

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

(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;

(4)

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

>	printlevel:=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;

(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);

(7)

>	Ntot:=ntot;

(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;

(9)

>

Nx;

(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));

(11)

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

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

>	Phi0:=Phi0+db:

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

(12)

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

(13)

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

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

(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;

(15)

>	phiaveAdd(Y00,N,M,Ny);

(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);

(17)

>

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

(18)

>	printlevel:=1:

>	Ymid:=copy(Y00):

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

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

(19)

>

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

(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;

(21)

>

V[nt];

(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;

(23)

>

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

>	tf:=TT[nt];

(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..NNh+h/2,-h/2..MMh+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]);

>

Download Phasefield2DBaseCodeParametric.mw

Calling UMFPACK sparsedirect linearsolve...

Asked: Dr. Venkat Subramanian 351 Product: Maple

February 16 2022

0 2

@acer,
In the past, you had shown how to use take a given sparse A matrix and store the LU components and use the same factorization for multiple backsolves at https://www.mapleprimes.com/posts/41191-Solving-Sparse-Linear-Systems-In-Maple#comment200817

When PDEs are solved, we are calling the A matrix at every time step (with some discretizations in x) and the factorization is done and stored at every time step (this is time-consuming and memory-consuming). Is it possible to call UMFPACK using only the sparse storage and entries? The main routine seems to provide the option

https://people.sc.fsu.edu/~jburkardt/f77_src/umfpack/umfpack.html

This would mean that the pattern is found only once for the Jacobian at t= 0, then only the non-zero sparse entries (vector/row, not a matrix) are updated at every time step.

To be clear, what I am asking for is create a random sparse matrix (say 4x4 or 10x10)
create a b Vector.

Solve AZ =b with method = SparseDirect. (This should internally create and store R, CC, X as at https://people.sc.fsu.edu/~jburkardt/f77_src/umfpack/umfpack.html)

Store R, CC, X and update only X for a different A matrix with the same sparsity pattern and solve for Z using stored (R, CC), and updated X.

Non symmetric sparse solver - non-implem...

Asked: Dr. Venkat Subramanian 351 Product: Maple

May 06 2020

2 1

Hi All,

https://www.maplesoft.com/support/helpJP/Maple/view.aspx?path=NAG%2ff11dec

This link shows that there are 3 options rgmres, cgs and bicgstab. When I use iterative solver, it uses the cgs option.

LinearAlgebra:-LinearSolve(
A, b, method= SparseIterative)

LinearSolve: "copying right hand side to enable external call"
LinearSolve: "using method" SparseIterative
LinearSolve: "using method " SparseIterative
LinearSolve: "calling external function"
LinearSolve: "using CGS method"
LinearSolve: "preconditioning with incomplete LU factorization"
LinearSolve: "level of fill = " 0
LinearSolve: "using complete pivoting strategy"
LinearSolve: "dimension of workspaces for preconditioner = " 44
LinearSolve: "using infinity norm in stopping criteria"
LinearSolve: "setting maximum iterations to " 200
LinearSolve: "setting tolerance to " .10e-7
LinearSolve: "NAG" hw_f11zaf
LinearSolve: "NAG" hw_f11daf
LinearSolve: "NAG" hw_f11dcf
LinearSolve: "number of iterations" 1
LinearSolve: "residual computed last as" 3.33066907387547e-016

How can I force Maple to use BiCGSTAB?

thanks

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

Venkat Subramanian

351 Reputation

13 Badges

MaplePrimes Activity

These are questions asked by

Inconsistent behavior with dsolve/numeri...

Will anyone be able to speed up this cod...

Will anyone be able to speed up this cod...

Calling UMFPACK sparsedirect linearsolve...

Non symmetric sparse solver - non-implem...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again