Questions - MaplePrimes

loops in MapleFlow...

Asked by: charlie_fcl 95

April 24 2022

3 2

Hi,

I like to know if it is possible to use a for loops in MapleFlow. I've tried a couple of times but it didn't work. For example,

n:=6
a:=Matrix(n,n)

for l from 1 to n do

   for c from 1 to n do

     a[l,c]:=2^(l+c)

   end do

end do

a=

doesn't produce any result for a

Best regards,

Numerical Scheme for Integro-differential equation...

Asked by: Tamour_Zubair 80

April 23 2022

1 1

Hello;

Hope you are fine. Can i apply numerical scheme on maple for the following problem. This in integro-differential equation i think. Waiting for kind response.

Thanks

"Help!for this series."...

Asked by: Saha 70

April 23 2022

0 1

could you please help me ,the maple code for this given series.

>

(1)

>

$"U[k+1](x):=solve((k+1)*U[k+1](x)+(x*(-1)^((k-1)/(2)))/(k!)-x^(2)*((e)^(x))/(10){6/(k!)-sum((2^(k[1]))/(k[1!])(5*delta[k[]-k[1]](x)+(2^(k-k[1])*(-1)^((k-k[1])/(2)))/((k-k[1])!)+(2^(k-k[1]+1)*(-1)^((k-k[1]-1)/(2)))/((k-k[1])!)),k[1]=0..k)}-(cos^(2)(x)+sin^()(x))*((&PartialD;)^2)/((&PartialD;)^( )x^2) [U[k](x)]-(e)^(x)[sum(1/(k[1]!){1/(k-k[1])(sum(sum(1/(k[3]!)*U[k[2]-k[3]](x)*U[k-k[1]-k[2]-1](x)},k[1]=0..k),k[2]=0..k-k[1]-1),k[3]=0..k[2])),U[k+1](x)]; od;"$

Error, unable to match delimiters

$"U[k+1](x):=solve((k+1)*U[k+1](x)+(x*(-1)^((k-1)/2))/(k!)-x^2*((e)^x)/10{6/(k!)-sum((2^(k[1]))/(k[1!])(5*delta[k-k[1]](x)+(2^(k-k[1])*(-1)^((k-k[1])/2))/((k-k[1])!)+(2^(k-k[1]+1)*(-1)^((k-k[1]-1)/2))/((k-k[1])!)),k[1]=0..k)}-(cos^2(x)+sin(x))*((&PartialD;)^2)/(&PartialD;x^2) [U[k](x)]-(e)^x[sum(1/(k[1]!){1/(k-k[1])(sum(sum(1/(k[3]!)*U[k[2]-k[3]](x)*U[k-k[1]-k[2]-1](x)},k[1]=0..k),k[2]=0..k-k[1]-1),k[3]=0..k[2])),U[k+1](x)]; od;"$

>

Download Chapter_6-Example-6.5.4.mw

How to make animation of matrix plotting ?...

Asked by: PavelGR 5

April 23 2022

1 2

How to make the matrix construction by points animated?
I sketched this code - it doesn't work(((

a := Matrix([stoimostyM, izmeniya]);
plots[animate](plot, [a[n, m]], n = 1 .. 2970, m = 1 .. 2, numpoints = 50, frames = 100);
Error, bad index into Matrix

Belousov–Zhabotinsky equation ...

Asked by: Pepini 255

April 23 2022

2 2

Hi! Do you know maybe how to solve equation with Laplace operator in Maple like BZ equation?

restart;
a := 0.75;
rho := u(t) + v(t) + w(t);
ode := diff(u(t), t) = 10*Delta(u(t)) + u*(-a*v - rho + 1);
ode1 := diff(v(t), t) = 10*Delta(v(t)) + v*(-a*w - rho + 1);
ode2 := diff(w(t), t) = 10*Delta(w(t)) + w*(-a*u - rho + 1);

Edit: Sorry I guess it should be function of three variables so u,v,w depends on (x,y,z) not strictle from time

I am wondering how to animate something like this from BZ equation:

How can I generate the following list from L=[x,y,...

Asked by: Chorux 60

April 23 2022

1 5

Given a list L=[x,y,z], what is the best way to generate the following sequential expression?
L₁=[[x],[y],[z],[x,x],[x,y],[[x,z],[y,y],[y,z],[z,z],[x,x,x],[x,x,y],[x,x,z],[x,y,y],[x,y,z],[y,y,y],[y,y,z],[y,z,z],[z,z,z]].

With L₁in mind (a list of list of at most three variables), how can one generate a list of lists of any number of variable from L.

I am currently out of options, a response will be highly appreciated.

Thanks

ERROR: Flow does not evaluate to a module...

Asked by: rs4670 10

April 22 2022

2 1

I have Maple 2021 and Mapleflow 2021 installed. The Maple 2021 is working fine. However, MapleFlow gives the error as follows. ERROR: Flow does not evaluate to a module. My system OS is Windows 10 updates till April 2022

Will anyone be able to speed up this code even mor...

Asked by: Dr. Venkat Subramanian 386

April 22 2022

3 17

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

Issue to make code of rk-2 for system of ODE's...

Asked by: Tamour_Zubair 80

April 21 2022

2 6

Hello;

Hope you are fine. I am devolping code for rk-2 at N=8 and further comparing my results with bultin comands on rk-2 on maple. My results is not matching with maple buliton commands. I am confused and tried a lot but not able to find the mistake. I have attached code. May be i have mistake in the calculation of K_1 and K_2. Kindly guide me please. Thanks a lot in advance.

Question_(1)_Posted.mw

''How do I solve a differential equation in Maple'...

Asked by: Saha 70

April 21 2022

2 3

I can't solve this series,please could you help me ,the maple code for this.

restart;
U[0](x):=(e)^(x);
exp(x)
n:=7: for k from 0 to n do U[k+1](x):=solve((k+1)U[k+1](x)-((&PartialD;)^2)/((&PartialD;)^( )x^2) U[k](x)-2 U[k](x)); end do;
Error, (in solve) cannot solve expressions with diff(diff(U[0](x), x), x) for x

can u help my homework...

Asked by: south 50

April 21 2022

2 3

Let the curves y=2x^3-x^2-5x and 𝑦=−𝑥^2+3𝑥be given.
a) Plot both curves on the same 𝑥y −axis.
b) Shade the region (between the curves).
c) Calculate the area of the region enclosed by the curves given.

How i can find on maple...

Asked by: south 50

April 21 2022

3 5

What values of a and b make f(x)=x³+ax²+bx have

a. A local maximum at x=−1 and a local minimum at x=3?

b. A local minimum at x=4 and a point of inflection at x=1?

How to solve ode system with equilibrium points...

Asked by: Prakash J 100

April 21 2022

2 2

eq1 := diff(x(t), t) = x(K[1]-x(t))-p*x(t-tau[1])*y(t); eq2 := diff(y(t), t) = y(K[2]-x(t))-q*y(t-tau[2])*x(t);
(x*, y*):= (p*K[2]-k[1])/(p*q-1), (q*K[1]-k[2])/(p*q-1);

where, k_1, k_2, p, q, tau_1, tau_2 are positive constants

Errors (red boxes) in Maple Flow...

Asked by: mapleflowuser 10

April 21 2022

2 1

I get this error when typing a sigma symbol with some subscripts.

Is there anyway to get rid of it?

If mapleflow doesn't allow typing in this way, how can I type it in an unexcutable Math format like document mode in Maple?

Help! Symbolic optimization...

Asked by: danortega 10

April 21 2022

2 3

Hi everyone! I'd really appreciate if I could get pointed in the right direction as I am a brand new maple user.

So im trying to solve this constrainted optimization problem (See picture) using Maple symbollically. I believe I should have a closed form solution given I can substitute the one constraint into the objective function. Specifically closed form solutions for the three phi variables.

Can someone point me in the right direction as to how I should go about this? I've already taken first order conditions and tried to using the solve() function to no avail, realizing my sytem of equations weren't linear );.

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

MaplePrimes Questions

loops in MapleFlow...

Numerical Scheme for Integro-differential equation...

"Help!for this series."...

How to make animation of matrix plotting ?...

Belousov–Zhabotinsky equation ...

How can I generate the following list from L=[x,y,...

ERROR: Flow does not evaluate to a module...

Will anyone be able to speed up this code even mor...

Issue to make code of rk-2 for system of ODE's...

''How do I solve a differential equation in Maple'...

can u help my homework...

How i can find on maple...

How to solve ode system with equilibrium points...

Errors (red boxes) in Maple Flow...

Help! Symbolic optimization...

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