Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

Problem 5.

Definition. The q-binomial coefficient ---denoted by qbinomial(n,k,q) --is defined by the formula
      qbinomial(n,k,q)    =       
"Product(q^i - 1, i = n - k + 1 .. n)/Product(q^i - 1, i = 1 .. 

   k)"

qbinomial(n,k,q) is a polynomial of degree k*(n-k) in the variable q. These polynomials are also called Gaussian polynomials or Gaussian  coefficients. They occur in various places in mathematics. For example qbinomial(n,k,q) is the number of k-dimensional subspaces of an n-dimensional vector space over a field with q elements.  An example:  qbinomial(6,3,q)  is the polynomial 

"q^9 + q^8 + 2*q^7 + 3*q^6 + 3*q^5 + 3*q^4 + 3*q^3 + 2*q^2 + q 

   + 1"

(a) Write a procedure to find qbinomial(n,k,q) and check to see that it agrees with the above for n=6  and k=3.  Note that to get nice output you will need to simplify, expand and sort the expression.

(b) Check whether or not the sequence of coefficients of the Gaussian polynomial qbinomial(n,k,q) is unimodal for  0 =< n < =10 and 0=< k < =n.  Do not print out all these in your final work. Just write a program to do the checking.  Have it print out the number of pairs (n,k) for which qbinomial(n,k,q) 's coefficients are not unimodal. This is a number from 0 to 66.

[Hint: You may use the following to convert the coefficients of a polynomial f in the variable q to a list:
>L:=[seq(coeff(f,q,i),i=0..degree(f,q))]; ]

Definition: A list or sequence of numbers 
                             "a[1]"


                             "a[2]"

, ..., 
                             "a[n]"

 is unimodal if there is an index i such that 
                         
                             "a[1]"

 <= 
                             "a[2]"

 <= . . . <= 
                             "a[i]"

 >= 
                           "a[i + 1]"

 >= 
                           "a[i + 2]"

 >= . . . >= 
                             "a[n]"


That is, the sequence is non-decreasing up to some point after which it is non-increasing. Note that i can be 1 or n.  A constant sequence is considered to be unimodal.

Examples of unimodal lists:

         [1, 1, 1, 1, 1],  
         [1,2,3,4,5,4,3,2,1], 
         [1, 2, 2, 3, 4, 5, 5, 5],   
         [5, 5, 4, 4, 3, 3, 1],   
         [1, 2, 2, 3, 3, 3, 4, 4, 2, 2, 1, 1, 1]

Examples of lists that are not unimodal:

          [1, 0, 1, 0], 
          [1, 1, 2, 2, 3, 4, 5, 2, 2, 6, 4, 2, 2, 1, 0]

(a) Write a procedure unimodal to check whether or not a list of numbers is unimodal. The input should be a list and the output should be true or false. 

One way to do this is to first write procedures called  increasing and decreasing which will check whether or not a sequence is non-decreasing (
                             "a[1]"

 <= 
                             "a[2]"

 <=...) or non-increasing (
                             "a[1]"

>= 
                             "a[2]"

>=...). Then for a list L use these procedures to check the two parts L[1..i] and L[i..n], n=nops(L) for each i. If you find an i such that the first list is "increasing" and the second is "decreasing" then you can return true. 

I have up to here.

I need help with this part.

For n from 10 to 15  check whether or not the sequence of binomial coefficients  

[binomial(n,0), binomial(n,1), binomial(n,2), . . ., binomial(n,n)]

is unimodal. 

Hello there, 

Would you allow me to ask this one question?

Is it possible to apply diffop2de() to three different functions in a single expression? For example, the attached worksheet shows an expression, 'eq_e5_9', where three time domain function derivatives are defined, D*Delta__delta(t), D*Delta__psi__d(t) and D*Delta__psi__q(t). However, once 'diffop2de(eq_e5_9, f(t))' is applied, none of them gets propely unrolled (e.g, D*Delta__delta(t)becomes d/dt*Delta__delta(t))

Besides, this is my third time attempting to post a question. The first one got disappeared with no sign, the second one came out with a site error message. 

In Kwon Park 

 

Download Q20220425.mw

Dear all

I have a nonlinear system of algebraic equations, I would like to solve it  without using fsolve and Newton's method. 
Maybe one can use, fixed point method or Broyden's method

Fixed_Broyden_method.mw

Can we solve the system of 4 equations, using the proposed methods

thank you 

Hi everyone, I have some questions about the printf func in maple.

As an image I attach below, I want to print this to the console so that I can presentation it for my teacher, because I'm using this palletes many times, only showing the result will make me and my teacher hard to check the right result. I trying using the sum palletes as a tring ("%s") and also the result as ("%f") too, but it wont work. Please help me with this, thanks a lot

Solve produces different output in the attachment depending on how it is used. Why is that and how can simplification to arctan(y/z) be avoided? Arctan(y/z) only gives correct angles for positive y and z.  I prefer arctan(y,z) output that I can subsequently simplify to the y and z ranges of interest (if possible). Imagine “wrong” simplification of complex algebraic output (e.g., from inverse kinematics).

Arctan.mw

I am using the method of alias(seq(c[k] = _C||k, k = 1..10)); for better latex of constants generated from solving an ode as recommended. See this for example. 

This works well 99.99% of the time. But now I noticed this in Maple 2022. Is this a display issue? When the constant is inside an inert Int it does not display the same as the other constant outside. Also the Latex is not the same. Even though lprint shows they are both correct.

This is the worksheet itself


 

restart;

interface(version);

`Standard Worksheet Interface, Maple 2022.0, Windows 10, March 8 2022 Build ID 1599809`

alias(seq(c[k] = _C||k, k = 1..10));
ode:=x*diff(y(x),x$2)-cos(x)*diff(y(x),x)+sin(x)*y(x)=2;
sol:=dsolve(ode);
lprint(sol)

c[1], c[2], c[3], c[4], c[5], c[6], c[7], c[8], c[9], c[10]

x*(diff(diff(y(x), x), x))-cos(x)*(diff(y(x), x))+sin(x)*y(x) = 2

y(x) = (c[2]+Int((c[1]+2*x)/(exp(Ci(x))*x^2), x))*exp(Ci(x))*x

y(x) = (c[2]+Int((c[1]+2*x)/exp(Ci(x))/x^2,x))*exp(Ci(x))*x

latex(sol)

y \! \left(x \right) =
\left(c_{2}+\textcolor{gray}{\int}\frac{\mathit{c[1]} +2 x}{{\mathrm e}^{\mathrm{Ci}\left(x \right)} x^{2}}\textcolor{gray}{d}x \right) {\mathrm e}^{\mathrm{Ci}\left(x \right)} x

restart;

ode:=x*diff(y(x),x$2)-cos(x)*diff(y(x),x)+sin(x)*y(x)=2;
sol:=dsolve(ode);
lprint(sol)

x*(diff(diff(y(x), x), x))-cos(x)*(diff(y(x), x))+sin(x)*y(x) = 2

y(x) = (_C2+Int((_C1+2*x)/(exp(Ci(x))*x^2), x))*exp(Ci(x))*x

y(x) = (_C2+Int((_C1+2*x)/exp(Ci(x))/x^2,x))*exp(Ci(x))*x

 


 

Download april_25_2022.mw

 

 

Hi,

I'm experimenting with some plot functions in Maple and cannot make the animate function to work, even for the help example (attached). It works fine when I insert it here but not in my worksheet.

with(plots)

NULL

plots[animate](plot, [A*sin(x), x = 0 .. 10], A = 0 .. 2)

 

NULL

Download animate_prob.mw

Any ideas?

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

 


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

restart

U[0](x) := x;

x

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

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:

Given a list L=[x,y,z], what is the best way to generate the following sequential expression?
L1=[[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 Lin 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

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

 

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


 

 

First 195 196 197 198 199 200 201 Last Page 197 of 2097