- Replies - MaplePrimes

Please try and let me know...

Answered: Dr. Venkat Subramanian 431

December 07 2022

0 0

If the solver works or fails for any index-1 DAE that you are testing.

Tridiagonal is easy...

Answered: Dr. Venkat Subramanian 431

August 27 2022

0 0

@Carl Love

The solver for the tridiagonal system is not relevant for general sparse systems that might arise for 2D models (block-banded). The only way to do that will be the alternating direction method which solves for multiple 1D models alternatively in x and y.

Read again...

Answered: Dr. Venkat Subramanian 431

August 25 2022

0 0

@Alger

Maple's Linearsolve was used

It is storing only the non-zero values of your matrix in vector format.

not an answer but see this thread...

Answered: Dr. Venkat Subramanian 431

August 25 2022

0 0

https://www.mapleprimes.com/questions/234050-Will-Anyone-Be-Able-To-Speed-Up-This

Typically, depending on the pattern, Pardiso can help you solve 1000000x1000000 very fast, but that is not directly shipped with Maple.
For umfpack if you have to solve only once, you can still do it, but you have to use the vector format (not the matrix format), see the discussion. It may be very involved depending on your system.

I believe we are done!...

Answered: Dr. Venkat Subramanian 431

May 31 2022

0 0

@ acer

Releasing the RAM using extfun and gc() and restarting Maple seems to work for reducing the RAM spike in the task manager.

Maple 2022 does this better than Maple 2020. This code works for N = 320 as well. The speed is comparable to PARDISO (slightly slower). Now, this code can be used by any Maple user! Thanks, acer.

>

restart:

>	Digits:=15:

>	with(LinearAlgebra):

>	t11:=time():

>	t12:=time[real]():

Code was updated on 05/28/22 to call the sparse jacobian in vector format. This enables compiled call for both the residues and the sparse entries of the Jacobian. The structure and 'sym' part of the UMFPACK linearsolver are called only once.

Code updated on 05/02/22 to minmize reduce overhead calls for LinearSolve (based on inputs from acer).
It directly calls MatVecSolve in UMFPACK. Many procedures are autocompiled.

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;MM:=NN;delta:=0.1;vel:=1.0;ki0:=1.0;tf:=2.0;

>

gc();

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

>

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:=evalf(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:

>

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

>

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.

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

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

>

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]),N::integer,M::integer)
local i::integer,j::integer,i1::integer,h::float[8];
option optimize, autocompile;
h:=1.0/(N-2):
for j from 1 to M do
for i from 1 to N do
i1:=i+(j-1)*N:
if i>1 and i <N and j>1 and j<M then
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])
-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):
end:
if i =1 and j =1 then ff[i1]:=-y[i1]:end:
if i =N and j =1 then ff[i1]:=-y[i1]:end:
if i =1 and j =M then ff[i1]:=-y[i1]:end:
if i =N and j =M then ff[i1]:=-y[i1]:end:

if i >1 and i<N and j =1 then ff[i1]:=-y[i1]+y[i1+N]:end:
if i >1 and i<N and j =M then ff[i1]:=-y[i1]+y[i1-N]+delta*h:end:
if j >1 and j<M and i =1 then ff[i1]:=-y[i1]+y[i1+1]:end:
if j >1 and j<M and i =N then ff[i1]:=-y[i1]+y[i1-1]:end:

od:
od:
end proc:

>

jsproc:=proc(Y00::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,jss::Vector(datatype=float[8]),N::integer,M::integer)
local count::integer,i::integer,j::integer,cent::float[8],top::float[8],bot::float[8],left::float[8], right::float[8],h::float[8];
option optimize, autocompile;
h:=1.0/(N-2):
count:=0:for j from 1 to M do
for i from 1 to N do
if i>1 and i<N and j>1 and j<M then
count:=count+1:
if j = 2 then bot:=1.0: else bot:=Y00[i,j]+Y00[i,j-1]:end:
jss[count]:=bot:

count:=count+1:
if i = 2 then left:=1.0; else left:=Y00[i,j]+Y00[i-1,j]:end:
jss[count]:=left:

count:=count+1:
cent:=-4*Y00[i,j]-Y00[i,j+1]-Y00[i,j-1]-Y00[i+1,j]
-Y00[i-1,j]-ki0*h*(1e-24+(Y00[i+1,j]-Y00[i-1,j])^2
+(Y00[i,j+1]-Y00[i,j-1])^2)^(1/2):
jss[count]:=cent:#print(cent);

count:=count+1:
if i = N-1 then
right:=1.0; else right:=Y00[i,j]+Y00[i+1,j]:end:
jss[count]:=right:

count:=count+1:
if j=M-1 then top:=1.0: else top:=Y00[i,j]+Y00[i,j+1]:end:
jss[count]:=top:
end:

if i =1 and j =1 then count:=count+1: cent:=-1.0: jss[count]:=cent: end:
if i =N and j =1 then count:=count+1: cent:=-1.0: jss[count]:=cent: end:
if i =1 and j =M then count:=count+1: cent:=-1.0: jss[count]:=cent: end:
if i =N and j =M then count:=count+1: cent:=-1.0: jss[count]:=cent: end:

if i >1 and i<N and j =1 then
count:=count+1:cent:=-1.0: jss[count]:=cent:
count:=count+1:top:=Y00[i,j]+Y00[i,j+1]:jss[count]:=top: end:

if i >1 and i<N and j =M then
count:=count+1:bot:=Y00[i,j]+Y00[i,j-1]: jss[count]:=bot:
count:=count+1:cent:=-1.0: jss[count]:=cent: end:

if j >1 and j<M and i =1 then
count:=count+1:cent:=-1.0: jss[count]:=cent:
count:=count+1:right:=Y00[i+1,j]+Y00[i,j]: jss[count]:=right:end:

if j >1 and j<M and i =N then
count:=count+1:left:=Y00[i,j]+Y00[i-1,j]: jss[count]:=left:
count:=count+1:cent:=-1.0: jss[count]:=cent:end:

od:
od:
end proc:

>

ApAiproc:=proc(Api::Vector(datatype=integer[4]),Aii::Vector(datatype=integer[4]),N::integer,M::integer)
local count::integer[8],i::integer,j::integer,ii::integer:
option optimize, autocompile;
count:=0:
for j from 1 to M do
for i from 1 to N do
ii:=i+(j-1)*N:
if i>1 and i<N and j>1 and j<M then
count:=count+1:
Aii[count]:=i+(j-1)*N-1-N:Api[ii]:=count-1:

count:=count+1:
Aii[count]:=i+(j-1)*N-1-1:

count:=count+1:
Aii[count]:=i+(j-1)*N-1:#print(cent);

count:=count+1:
Aii[count]:=i+(j-1)*N-1+1:

count:=count+1:
Aii[count]:=i+(j-1)*N-1+N:
end:

if i =1 and j =1 then count:=count+1: Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:end:
if i =N and j =1 then count:=count+1: Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:end:
if i =1 and j =M then count:=count+1: Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:end:
if i =N and j =M then count:=count+1: Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:end:

if i >1 and i<N and j =1 then
count:=count+1:Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:
count:=count+1:Aii[count]:=i+(j-1)*N-1+N:end:

if i >1 and i<N and j =M then
count:=count+1:Aii[count]:=i+(j-1)*N-1-N:Api[ii]:=count-1:
count:=count+1:Aii[count]:=i+(j-1)*N-1:end:

if j >1 and j<M and i =1 then
count:=count+1:Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:
count:=count+1:Aii[count]:=i+(j-1)*N-1+1:end:

if j >1 and j<M and i =N then
count:=count+1:Aii[count]:=i+(j-1)*N-1-1:Api[ii]:=count-1:
count:=count+1:Aii[count]:=i+(j-1)*N-1:end:

od:
od:
ii:=ii+1:
Api[ii]:=Api[ii-1]+1:
end proc:

>	#Compiler:-Compile(ApAiproc);

The Jacobian for the residues is coded in the procedure Jac1.Note that in this code the full matrix for the Jacobian is not used and Jac1 is not used. Jac1 is provided only to help the readers understand the procedure jsproc (nonzero values of Jac1) and the sparse storage format

>

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

for j from 2 to M-1 do
for i from 2 to N-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:
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:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*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:
od:od:
NULL;
end proc:

>

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];
option optimize, autocompile;
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];
option optimize, autocompile;
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];
option optimize, autocompile;
e1:=1e-6:
nx:=0.0:
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):
F0[i,j]:=nx:
od:od:
end proc:

>

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

>	phiaveAdd:=proc(Y00::Matrix(datatype=float[8]),N::integer,M::integer,Ny::integer) local i::integer,j::integer,phiave::float; option optimize, autocompile; 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:

>

EFAdd:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>

EFAdd2:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>

EFAdd3:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>	jtot:=(N-2)42+4+(N-2)(M-2)5:

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

>	#j00:=Matrix(1..ntot,1..ntot,datatype=float[8],storage=sparse): F0 :=Matrix(1..N,1..M,datatype=float[8]): Phi0:=Vector(1..ntot,datatype=float[8]): ff:=copy(Phi0):

>	Aii:=Vector(jtot,datatype=integer[4]):Api:=Vector(ntot+1,datatype=integer[4]): jss:=Vector(jtot,datatype=float[8]):

>	ApAiproc(Api,Aii,N,M);

>	Api:=convert(Api,Vector,datatype=integer[8]):Aii:=convert(Aii,Vector,datatype=integer[8]):

>	evalf(Eqs11(Y00,Phi0,delta,ki0,ff,N,M)):

>	jsproc(Y00,Phi0,delta,ki0,jss,N,M);

>	#Jac1(Y00,Phi0,delta,ki0,j00,N,M):

>	#with(LinearAlgebra):

>	umfsym:=define_external("umfpack_dl_symbolic", 'm'::(integer[8]), 'n'::(integer[8]), 'Ap'::ARRAY(integer[8]), 'Ai'::ARRAY(integer[8]), 'Ax'::ARRAY(float[8]), 'Sym'::REF(integer[8]), 'Control'::ARRAY(float[8]), 'Info'::ARRAY(float[8]), LIB="umfpack_dlong"):

>

>	umfnum:=define_external("umfpack_dl_numeric", 'Ap'::ARRAY(integer[8]), 'Ai'::ARRAY(integer[8]), 'Ax'::ARRAY(float[8]), 'Sym'::(integer[8]), 'Num'::REF(integer[8]), 'Control'::ARRAY(float[8]), 'Info'::ARRAY(float[8]), LIB="umfpack_dlong"):

>	umfdefaults:=define_external("umfpack_dl_defaults", 'Control'::ARRAY(float[8]), LIB="umfpack_dlong"):

>	Control:=Vector[row](1..20,datatype=float[8]):

>	umfdefaults(Control);

>	#umfdefaults;

>	seq(Control[i],i=1..20);

HFloat(1.0), HFloat(0.2), HFloat(0.2), HFloat(0.1), HFloat(32.0), HFloat(0.0), HFloat(0.7), HFloat(2.0), HFloat(0.0), HFloat(0.0), HFloat(0.0), HFloat(0.0), HFloat(0.01), HFloat(0.0), HFloat(10.0), HFloat(0.0010), HFloat(1.0), HFloat(0.5), HFloat(0.0), HFloat(1.0)

>	#Control[5]:=64.;

>	#seq(Control[i],1=1..4);

>	Info:=Vector[row](1..91,datatype=float[8]):

>	extlib:=ExternalCalling:-ExternalLibraryName("linalg",'HWFloat'): umfSolve:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_MatVecSolve',extlib):

>	umfsym(ntot,ntot, Api, Aii, jss, 'Sym', Control, Info):

>	umfnum(Api, Aii, jss, Sym, 'Num', Control, Info):

>	db:=Copy(Phi0):

>	db:=umfSolve(ntot,jtot,Num,-ff):

>	extfun:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_FreeNumeric', extlib):

>	PhiAdd(ntot,Phi0,db):

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

>	#extfun(Num);

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

>

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

>

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

>	YdatStore:=proc(Y00::Matrix(datatype=float[8]),Ydat::Array(datatype=float[8]),N::integer,M::integer,jj::integer) local i::integer,j::integer; option optimize, autocompile; 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:

>

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

>	Ymid:=copy(Y00):

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

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

>

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

>	HF:=eval(WENO3):

>

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

>

ii:=0:
while tt<tf do
HF(Y00,Phi0,F0,evalf(dt),N,M,delta):
EFAdd(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
Eqs11(Ymid,Phi0,delta,ki0,ff,N,M);
jsproc(Ymid,Phi0,delta,ki0,jss,N,M);
umfnum(Api, Aii, jss, Sym, 'Num', Control, Info);
db:=umfSolve(ntot,jtot,Num,-ff):extfun('Num');Num:=-1:
PhiAdd(ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd2(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
Eqs11(Ymid,Phi0,delta,ki0,ff,N,M);
jsproc(Ymid,Phi0,delta,ki0,jss,N,M);
umfnum(Api, Aii, jss, Sym, 'Num', Control, Info);
db:=umfSolve(ntot,jtot,Num,-ff):extfun('Num');Num:=-1:
PhiAdd(ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd3(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
Eqs11(Y00,Phi0,delta,ki0,ff,N,M);
jsproc(Y00,Phi0,delta,ki0,jss,N,M);
umfnum(Api, Aii, jss, Sym, 'Num', Control, Info);
db:=umfSolve(ntot,jtot,Num,-ff):extfun('Num');Num:=-1:
PhiAdd(ntot,Phi0,db):
ii:=ii+1:
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:
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);
#YdatStore(Y00,Ydat,N,M,ii+1);
Phiave[ii]:=phiaveAdd(Y00,N,M,MM);
#gc();
end:

>

nt:=ii;

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,axes=boxed);

>	V[0];V[2];

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

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

>

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

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,axes=boxed);

>

>	save Y00,"Y0data.m";

>	save tf,"tdata.m";

>

Download Phasefield2DBaseCodeParametricUMFPACKoptimizeVS.mw

Faster code almost done with UMFPACK so ...

Answered: Dr. Venkat Subramanian 431

May 30 2022

0 0

@acer

I was able to write the analytic jacobian in sparse storage vector format and compile it. This makes the code run faster and is probably comparable to PARDISO at lower N values (100 or so). For N = 100, it takes 4 to 5 seconds. But at higher N values, UMFPACK becomes slower with PARDISO doing a better job using all the cores in the computer.
For some weird reason, this approach is showing higher RAM usage in the task manager. Since LinearSolve based on Jacobian Matrix call was slower, I was not able to check its RAM use. But it is possible that the Control option is better handled by the inbuit LinearSolve for better memory usage. I am not sure. Note that my code is not showing high memory usage in the Maple window, but Java seems to be eating up RAM in the task manager.

>

restart:

>	Digits:=15:

>	with(LinearAlgebra):

>	t11:=time():

>	t12:=time[real]():

Code was updated on 05/28/22 to call the sparse jacobian in vector format. This enables compiled call for both the residues and the sparse entries of the Jacobian. The structure and 'sym' part of the UMFPACK linearsolver are called only once.

Code updated on 05/02/22 to minmize reduce overhead calls for LinearSolve (based on inputs from acer).
It directly calls MatVecSolve in UMFPACK. Many procedures are autocompiled.

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;MM:=NN;delta:=0.1;vel:=1.0;ki0:=1.0;tf:=2.0;

>

gc();

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

>

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:=evalf(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:

>

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

>

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.

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

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

>

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]),N::integer,M::integer)
local i::integer,j::integer,i1::integer,h::float[8];
option optimize, autocompile;
h:=1.0/(N-2):
for j from 1 to M do
for i from 1 to N do
i1:=i+(j-1)*N:
if i>1 and i <N and j>1 and j<M then
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])
-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):
end:
if i =1 and j =1 then ff[i1]:=-y[i1]:end:
if i =N and j =1 then ff[i1]:=-y[i1]:end:
if i =1 and j =M then ff[i1]:=-y[i1]:end:
if i =N and j =M then ff[i1]:=-y[i1]:end:

if i >1 and i<N and j =1 then ff[i1]:=-y[i1]+y[i1+N]:end:
if i >1 and i<N and j =M then ff[i1]:=-y[i1]+y[i1-N]+delta*h:end:
if j >1 and j<M and i =1 then ff[i1]:=-y[i1]+y[i1+1]:end:
if j >1 and j<M and i =N then ff[i1]:=-y[i1]+y[i1-1]:end:

od:
od:
end proc:

>

jsproc:=proc(Y00::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,jss::Vector(datatype=float[8]),N::integer,M::integer)
local count::integer,i::integer,j::integer,cent::float[8],top::float[8],bot::float[8],left::float[8], right::float[8],h::float[8];
option optimize, autocompile;
h:=1.0/(N-2):
count:=0:for j from 1 to M do
for i from 1 to N do
if i>1 and i<N and j>1 and j<M then
count:=count+1:
if j = 2 then bot:=1.0: else bot:=Y00[i,j]+Y00[i,j-1]:end:
jss[count]:=bot:

count:=count+1:
if i = 2 then left:=1.0; else left:=Y00[i,j]+Y00[i-1,j]:end:
jss[count]:=left:

count:=count+1:
cent:=-4*Y00[i,j]-Y00[i,j+1]-Y00[i,j-1]-Y00[i+1,j]
-Y00[i-1,j]-ki0*h*(1e-24+(Y00[i+1,j]-Y00[i-1,j])^2
+(Y00[i,j+1]-Y00[i,j-1])^2)^(1/2):
jss[count]:=cent:#print(cent);

count:=count+1:
if i = N-1 then
right:=1.0; else right:=Y00[i,j]+Y00[i+1,j]:end:
jss[count]:=right:

count:=count+1:
if j=M-1 then top:=1.0: else top:=Y00[i,j]+Y00[i,j+1]:end:
jss[count]:=top:
end:

if i =1 and j =1 then count:=count+1: cent:=-1.0: jss[count]:=cent: end:
if i =N and j =1 then count:=count+1: cent:=-1.0: jss[count]:=cent: end:
if i =1 and j =M then count:=count+1: cent:=-1.0: jss[count]:=cent: end:
if i =N and j =M then count:=count+1: cent:=-1.0: jss[count]:=cent: end:

if i >1 and i<N and j =1 then
count:=count+1:cent:=-1.0: jss[count]:=cent:
count:=count+1:top:=Y00[i,j]+Y00[i,j+1]:jss[count]:=top: end:

if i >1 and i<N and j =M then
count:=count+1:bot:=Y00[i,j]+Y00[i,j-1]: jss[count]:=bot:
count:=count+1:cent:=-1.0: jss[count]:=cent: end:

if j >1 and j<M and i =1 then
count:=count+1:cent:=-1.0: jss[count]:=cent:
count:=count+1:right:=Y00[i+1,j]+Y00[i,j]: jss[count]:=right:end:

if j >1 and j<M and i =N then
count:=count+1:left:=Y00[i,j]+Y00[i-1,j]: jss[count]:=left:
count:=count+1:cent:=-1.0: jss[count]:=cent:end:

od:
od:
end proc:

>

ApAiproc:=proc(Api::Vector(datatype=integer[4]),Aii::Vector(datatype=integer[4]),N::integer,M::integer)
local count::integer[8],i::integer,j::integer,ii::integer:
option optimize, autocompile;
count:=0:
for j from 1 to M do
for i from 1 to N do
ii:=i+(j-1)*N:
if i>1 and i<N and j>1 and j<M then
count:=count+1:
Aii[count]:=i+(j-1)*N-1-N:Api[ii]:=count-1:

count:=count+1:
Aii[count]:=i+(j-1)*N-1-1:

count:=count+1:
Aii[count]:=i+(j-1)*N-1:#print(cent);

count:=count+1:
Aii[count]:=i+(j-1)*N-1+1:

count:=count+1:
Aii[count]:=i+(j-1)*N-1+N:
end:

if i =1 and j =1 then count:=count+1: Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:end:
if i =N and j =1 then count:=count+1: Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:end:
if i =1 and j =M then count:=count+1: Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:end:
if i =N and j =M then count:=count+1: Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:end:

if i >1 and i<N and j =1 then
count:=count+1:Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:
count:=count+1:Aii[count]:=i+(j-1)*N-1+N:end:

if i >1 and i<N and j =M then
count:=count+1:Aii[count]:=i+(j-1)*N-1-N:Api[ii]:=count-1:
count:=count+1:Aii[count]:=i+(j-1)*N-1:end:

if j >1 and j<M and i =1 then
count:=count+1:Aii[count]:=i+(j-1)*N-1:Api[ii]:=count-1:
count:=count+1:Aii[count]:=i+(j-1)*N-1+1:end:

if j >1 and j<M and i =N then
count:=count+1:Aii[count]:=i+(j-1)*N-1-1:Api[ii]:=count-1:
count:=count+1:Aii[count]:=i+(j-1)*N-1:end:

od:
od:
ii:=ii+1:
Api[ii]:=Api[ii-1]+1:
end proc:

>	#Compiler:-Compile(ApAiproc);

The Jacobian for the residues is coded in the procedure Jac1.Note that in this code the full matrix for the Jacobian is not used and Jac1 is not used. Jac1 is provided only to help the readers understand the procedure jsproc (nonzero values of Jac1) and the sparse storage format

>

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

for j from 2 to M-1 do
for i from 2 to N-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:
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:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*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:
od:od:
NULL;
end proc:

>

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];
option optimize, autocompile;
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];
option optimize, autocompile;
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];
option optimize, autocompile;
e1:=1e-6:
nx:=0.0:
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):
F0[i,j]:=nx:
od:od:
end proc:

>

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

>	phiaveAdd:=proc(Y00::Matrix(datatype=float[8]),N::integer,M::integer,Ny::integer) local i::integer,j::integer,phiave::float; option optimize, autocompile; 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:

>

EFAdd:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>

EFAdd2:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>

EFAdd3:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>	jtot:=(N-2)42+4+(N-2)(M-2)5:

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

>	#j00:=Matrix(1..ntot,1..ntot,datatype=float[8],storage=sparse): F0 :=Matrix(1..N,1..M,datatype=float[8]): Phi0:=Vector(1..ntot,datatype=float[8]): ff:=copy(Phi0):

>	Aii:=Vector(jtot,datatype=integer[4]):Api:=Vector(ntot+1,datatype=integer[4]): jss:=Vector(jtot,datatype=float[8]):

>	ApAiproc(Api,Aii,N,M);

>	Api:=convert(Api,Vector,datatype=integer[8]):Aii:=convert(Aii,Vector,datatype=integer[8]):

>	evalf(Eqs11(Y00,Phi0,delta,ki0,ff,N,M)):

>	jsproc(Y00,Phi0,delta,ki0,jss,N,M);

>	#Jac1(Y00,Phi0,delta,ki0,j00,N,M):

>	#with(LinearAlgebra):

>	umfsym:=define_external("umfpack_dl_symbolic", 'm'::(integer[8]), 'n'::(integer[8]), 'Ap'::ARRAY(integer[8]), 'Ai'::ARRAY(integer[8]), 'Ax'::ARRAY(float[8]), 'Sym'::REF(integer[8]), 'Control'::ARRAY(float[8]), 'Info'::ARRAY(float[8]), LIB="umfpack_dlong"):

>

>	umfnum:=define_external("umfpack_dl_numeric", 'Ap'::ARRAY(integer[8]), 'Ai'::ARRAY(integer[8]), 'Ax'::ARRAY(float[8]), 'Sym'::(integer[8]), 'Num'::REF(integer[8]), 'Control'::ARRAY(float[8]), 'Info'::ARRAY(float[8]), LIB="umfpack_dlong"):

>	umfdefaults:=define_external("umfpack_dl_defaults", 'Control'::ARRAY(float[8]), LIB="umfpack_dlong"):

>	Control:=Vector[row](1..20,datatype=float[8]):

>	umfdefaults(Control);

>	#umfdefaults;

>	seq(Control[i],i=1..20);

>	#Control[5]:=64.;

>	#seq(Control[i],1=1..4);

>	Info:=Vector[row](1..91,datatype=float[8]):

>	extlib:=ExternalCalling:-ExternalLibraryName("linalg",'HWFloat'): umfSolve:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_MatVecSolve',extlib):

>	umfsym(ntot,ntot, Api, Aii, jss, 'Sym', Control, Info):

>	umfnum(Api, Aii, jss, Sym, 'Num', Control, Info):

>	db:=Copy(Phi0):

>	db:=umfSolve(ntot,jtot,Num,-ff):

>	#extfun:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_FreeNumeric', #extlib):

>	PhiAdd(ntot,Phi0,db):

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

>	#extfun(Num);

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

>

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

>

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

>	YdatStore:=proc(Y00::Matrix(datatype=float[8]),Ydat::Array(datatype=float[8]),N::integer,M::integer,jj::integer) local i::integer,j::integer; option optimize, autocompile; 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:

>

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

>	Ymid:=copy(Y00):

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

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

>

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

>	HF:=eval(WENO3):

>

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

>

ii:=0:
while tt<tf do
HF(Y00,Phi0,F0,evalf(dt),N,M,delta):
EFAdd(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
Eqs11(Ymid,Phi0,delta,ki0,ff,N,M);
jsproc(Ymid,Phi0,delta,ki0,jss,N,M);
umfnum(Api, Aii, jss, Sym, 'Num', Control, Info);
db:=umfSolve(ntot,jtot,Num,-ff):#extfun(Num);Num:=-1:
PhiAdd(ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd2(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
Eqs11(Ymid,Phi0,delta,ki0,ff,N,M);
jsproc(Ymid,Phi0,delta,ki0,jss,N,M);
umfnum(Api, Aii, jss, Sym, 'Num', Control, Info);
db:=umfSolve(ntot,jtot,Num,-ff):#extfun(Num);Num:=-1:
PhiAdd(ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd3(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
Eqs11(Y00,Phi0,delta,ki0,ff,N,M);
jsproc(Y00,Phi0,delta,ki0,jss,N,M);
umfnum(Api, Aii, jss, Sym, 'Num', Control, Info);
db:=umfSolve(ntot,jtot,Num,-ff):#extfun(Num);Num:=-1:
PhiAdd(ntot,Phi0,db):
ii:=ii+1:
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:
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);
#YdatStore(Y00,Ydat,N,M,ii+1);
Phiave[ii]:=phiaveAdd(Y00,N,M,MM);
gc();
end:

>

nt:=ii;

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,axes=boxed);

>	V[0];V[2];

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

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

>

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

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,axes=boxed);

>

>	save Y00,"Y0data.m";

>	save tf,"tdata.m";

>

Download Phasefield2DBaseCodeParametricUMFPACKoptimizeVS.mw

fixed this...

Answered: Dr. Venkat Subramanian 431

May 06 2022

0 0

@ acer
Removing 'WRAPPER' made it work.

Updated code...

Answered: Dr. Venkat Subramanian 431

May 02 2022

0 0

>

restart:

>	Digits:=15:

>

>	t11:=time():

>	t12:=time[real]():

>

Code updated on 05/02/22, making use of MatVecSolve in UMFPACK. Many procedures are compiled using autocompile.
Code updated on 04/26/22 to minmize reduce overhead calls for LinearSolve (based on inputs from acer).

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;MM:=NN;delta:=0.1;vel:=1.0;ki0:=1.0;tf:=2.0;

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

>

SLU:=proc(A::Matrix(datatype=float[8]),x::Vector(datatype=float[8]))
local extlib,HWcall1,HWcall2,HWcall3,Anz,NumericA,nrowsA,ncolsA,res;
option optimize;
extlib:=ExternalCalling:-ExternalLibraryName("linalg",'HWFloat'):
HWcall1:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_MatFactor',extlib):
HWcall2:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_MatVecSolve',extlib):
HWcall3:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_FreeNumeric',extlib):
nrowsA,ncolsA:=op(1,A);#print(nrowsA);
Anz,NumericA := HWcall1(nrowsA,ncolsA,A);
res := HWcall2(nrowsA,Anz,NumericA,x);
HWcall3(NumericA);
res;
end proc:

>

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:=evalf(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:

>

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

>

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.

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

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

>

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]),N::integer,M::integer)
local i::integer,j::integer,i1::integer,h::float[8];
option optimize, autocompile;
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:
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],storage=sparse),N::integer,M::integer)
local i::integer,j::integer,i1::integer,h::float[8];
#option optimize, autocompile;
option optimize;
h:=1.0/(N-2):

for j from 2 to M-1 do
for i from 2 to N-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:
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:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*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:
od:od:
NULL;
end proc:

>

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];
option optimize, autocompile;
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];
option optimize, autocompile;
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];
option optimize, autocompile;
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):
F0[i,j]:=nx:
od:od:
end proc:

>

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

>	phiaveAdd:=proc(Y00::Matrix(datatype=float[8]),N::integer,M::integer,Ny::integer) local i::integer,j::integer,phiave::float; option optimize, autocompile; 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:

>

EFAdd:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>

EFAdd2:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>

EFAdd3:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,vel::float,ki0::float,N::integer,M::integer)
local i::integer,j::integer;
option optimize, autocompile;
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:

>

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

>	j00:=Matrix(1..ntot,1..ntot,datatype=float[8],storage=sparse): F0 :=Matrix(1..N,1..M,datatype=float[8]): Phi0:=Vector(1..ntot,datatype=float[8]): ff:=copy(Phi0):

>

>	evalf(Eqs11(Y00,Phi0,delta,ki0,ff,N,M)):

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

>	db:=SLU(j00,-ff):

>	PhiAdd(ntot,Phi0,db):

>

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

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

>

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

>

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

>	YdatStore:=proc(Y00::Matrix(datatype=float[8]),Ydat::Array(datatype=float[8]),N::integer,M::integer,jj::integer) local i::integer,j::integer; option optimize, autocompile; 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:

>

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

>	Ymid:=copy(Y00):

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

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

>

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.

>

ii:=0:
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,N,M));
Jac1(Ymid,Phi0,delta,ki0,j00,N,M):
db:=SLU(j00,-ff):
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,N,M));
Jac1(Ymid,Phi0,delta,ki0,j00,N,M):
db:=SLU(j00,-ff):
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,N,M));
Jac1(Y00,Phi0,delta,ki0,j00,N,M):
db:=SLU(j00,-ff):
PhiAdd(ntot,Phi0,db):
ii:=ii+1:
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:
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);
YdatStore(Y00,Ydat,N,M,ii+1);
Phiave[ii]:=phiaveAdd(Y00,N,M,MM);
end:

>

nt:=ii;

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

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

>

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,axes=boxed);

>

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

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,axes=boxed);

>

>	save Y00,"Y0data.m";

>	save tf,"tdata.m";

>

Download Phasefield2DBaseCodeParametricMay02.mw

@acer

My student TaeJin was able to modify the call to UMFPACK by calling hw_SpUMFPACK_MatVecSolve instead of hw_SpUMFPACK_MatMatSolve. This avoids the agony of converting vector to matrix to vector after linearsolver call. The code is slightly faster.

In addition, the use of autocompile makes the code cleaner, even though it might make the codes slightly slower sometimes.

thanks...

Answered: Dr. Venkat Subramanian 431

April 30 2022

0 0

@acer

looking forward to it.

bump up...

Answered: Dr. Venkat Subramanian 431

April 27 2022

0 0

@acer

bump up. Any luck with this? Is it possible to store the compiled procedures in a library or a code so that they can be used for future uses?

thanks

VS

of course Linux is ok to start with...

Answered: Dr. Venkat Subramanian 431

April 23 2022

0 0

@ acer

I am ok if you start with Linux. We can always create codes and run them later in linux. But we have to provide a way for the users to run the code with pardiso solver (and umfpack solver) in a smooth manner so that the codes can be used without too much frustration.

After 10^4*10^4 iterated solvers may be faster, but direct solvers are very robust, and may be needed for these problems.

Windows please...

Answered: Dr. Venkat Subramanian 431

April 23 2022

0 0

@acer

Thanks. I was able to run your code. The wrapper you wrote was very helpful for me earlier for Rosenbrock methods (LU decomposition is done once with multiple solves).
If you run pardiso codes for N =160, 200 or more, you will see more gain from pardiso compared to the umfpack code.

I ran the standard code (umfpack code) in Maple 2020 (it works in 2021 or 2022 as well).

The pardiso code was run in Maple 2021.1 in windows version.

I prefer the Windows version, as I start building my codes using the classic worksheet Maple (in Maple 14 or so) and then move to mw format.

For the SLU wrapper, why do we need SLU can't we directly use

Anz,NumericA := HWcall[1](nrowsA,ncolsA,A);
res := HWcall[2](nrowsA,m,Anz,NumericA,localx);

after defing HWcall[1], call[2], etc?

thanks a lot

VS

2 of the 3 suggestions...

Answered: Dr. Venkat Subramanian 431

April 22 2022

0 0

@ acer,

I am not sure. I believe I may have already implemented the first 2 of your 3 suggestions in the pardiso code.

Excellent...

Answered: Dr. Venkat Subramanian 431

April 22 2022

0 0

@acer

Awesome. Please do so. Thanks. Phase-field models might require 500 or more nodes in each dimension for sufficient precision. Please find enclosed a pardiso based linear solver code (in which I had optimized the calls to the sparse matrix storage using just the non-zero entries). This code will not directly run in Maple, unless you get the MKL specific dll setup properly, but you can see the structure of the code here. I strongly believe UMFPACK can be used similarly (based on the instructions given by UMFPACK developers).

In the paper, we also introduce new RK methods called approximate RK3a in which linearsolver is called only once during the first stage, with interpolations at the other stages.

The Pardiso based Maple code is attached below. We were able to run 1000x1000 or more in an hour or so. This code also defines all the necessary procedures and uses the same compiled procedures to simulate the results for a different number of node points.

Most of these procedures can be preloaded as a separate code (or module) and directly called with N and design parameters in the main Maple window. But I want to give the codes with a minimum number of commands and full transparency.

>

restart:

>

04/21/2022. Code written by Dr. Venkat Subramanian and MAPLE group at UT. The code is provided as open-access with no restrictions.

This code is very similar to the standard 2D Battphase code. The only difference is the Pardiso linear solver is called after changing the dll files. Both the clock time and total computation time are reported. CPU reported in Table 3 at the end of the code for 160 node points in each direction. The entire Jacobian is not used. Only the non-zero sparse entries are identified (only once) and updated with time.

>	t11:=time[real]():t12:=time():

>

gc();

>	Digits:=15;

(1)

>

>	interface(prettyprint=2);

(2)

>

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],N::integer,M::integer,e1::float[8])
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(e1,Y00[i,j]-dt*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],N::integer,M::integer,e1::float[8])
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(e1,Y00[i,j]*3/4.+Ymid[i,j]/4.-dt/4.*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],N::integer,M::integer,e1::float[8])
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(e1,Y00[i,j]*1/3.+Ymid[i,j]*2/3.-dt*2/3.*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):

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

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

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

>	PhiSwitch:=Compiler:-Compile(PhiSwitch):

>	MidPred:=proc(N::integer,Phi0::Vector(datatype=float[8]),Phitemp::Vector(datatype=float[8]),Phimid::Vector(datatype=float[8])) local i::integer; for i from 1 to N do Phimid[i]:=Phitemp[i]/2.0+Phi0[i]/2.0:od: end proc:

>	MidPred:=Compiler:-Compile(MidPred):

>	y0proc:=proc(NN::integer,Y00::Matrix(datatype=float[8]),beta::float[8],e1::float[8]) local i,j,xx,yy,N,h,w,MM,M,Ny,rf,f,ff,rr;

>	N:=NN+2:h:=1.0/NN:w:=h*beta:

>	MM:=NN;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;
#if j<4 and i <N/5+1 then Y00[i,j]:=1e-9 else Y00[i,j]:=1.0;end:
#Y00[i,j]:=1.0:
#if rr <9/100. then Y00[i,j]:=1e-12; else Y00[i,j]:=1.0:end:
#if rr=9/100. then Y00[i,j]:=0.5:end:
Y00[i,j]:=max(e1,0.5+0.5*tanh((sqrt(rr)-0.3)/w/sqrt(2.0))):
#if abs(xx-0.5)>=0.25 and yy<=0.5 then Y00[i,j]:=1e-9 else Y00[i,j]:=1.0:end:
#if yy<=0.1 then Y00[i,j]:=1e-9 :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:

>

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

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

>	MM:=NN;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;
#if j<4 and i <N/5+1 then Y00[i,j]:=1e-9 else Y00[i,j]:=1.0;end:
#Y00[i,j]:=1.0:
#if rr <9/100. then Y00[i,j]:=1e-12; else Y00[i,j]:=1.0:end:
#if rr=9/100. then Y00[i,j]:=0.5:end:
#Y00[i,j]:=max(1e-12,0.5+0.5*tanh((sqrt(rr)-0.3)/w/sqrt(2.0))):
if abs(xx-0.5)>=0.25 and yy<=0.5 then Y00[i,j]:=1e-9 else Y00[i,j]:=1.0:end:
if yy<=0.1 then Y00[i,j]:=1e-9 :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:

>	y0proc:=Compiler:-Compile(y0proc):

>	printlevel:=2;

(3)

>

Eqs11:=proc(N::integer,M::integer,Y0::Matrix(datatype=float[8]),F::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),v0::float,ff::Vector(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
h:=1.0/(N-2):
j:=1: for i from 1 to N do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+N]:
od:

for j from 2 to M-1 do
i:=1:
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+1]:
for i from 2 to N-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])
-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:
i:=N:
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-1]+0*h:
od:

j:=M:
for i from 1 to N do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-N]+v0*h:
od:

end proc:

>	Eqs11:=Compiler:-Compile(Eqs11):

>

XX:=proc(N::integer,M::integer,Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),v0::float,x0::Vector(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8],count::integer;
h:=1.0/(N-2):
count:=0:
j:=1:
for i from 1 to N do
i1:=i+(j-1)*N:
#ff[i1]:=-y[i1]+y[i1+N]:
count:=count+1:
x0[count]:=1.0:
count:=count+1:
x0[count]:=-1.0:
#j00[i1,i1]:=1.0:j00[i1,i1+N]:=-1.00:
od:

for j from 2 to M-1 do
i:=1:
i1:=i+(j-1)*N:
count:=count+1:
x0[count]:=1.0:
count:=count+1:
x0[count]:=-1.0:
#j00[i1,i1]:=1.0:j00[i1,i1+1]:=-1.00:
for i from 2 to N-1 do
i1:=i+(j-1)*N:
count:=count+1:
x0[count]:=-Y0[i,j-1]-Y0[i,j]:
count:=count+1:
x0[count]:=-Y0[i-1,j]-Y0[i,j]:
count:=count+1:
x0[count]:=4*Y0[i,j]+Y0[i,j+1]+Y0[i,j-1]+Y0[i+1,j]+Y0[i-1,j]+h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
count:=count+1:
x0[count]:=-Y0[i+1,j]-Y0[i,j]:
count:=count+1:
x0[count]:=-Y0[i,j+1]-Y0[i,j]:
#j00[i1,i1]:=4*Y0[i,j]+Y0[i,j+1]+Y0[i,j-1]+Y0[i+1,j]+Y0[i-1,j]+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:
i:=N:
count:=count+1:
x0[count]:=-1.0:
count:=count+1:
x0[count]:=1.0:
i1:=i+(j-1)*N:
#j00[i1,i1]:=1.0:j00[i1,i1-1]:=-1.00:
od:

j:=M:
for i from 1 to N do
i1:=i+(j-1)*N:
count:=count+1:
x0[count]:=-1.0:
count:=count+1:
x0[count]:=1.0:
#j00[i1,i1]:=1.0:j00[i1,i1-N]:=-1.00:
#ff[i1]:=-y[i1]+y[i1-N]+v0*h:
od:

end proc:

>	XX:=Compiler:-Compile(XX):

>	CBproc:=proc(N::integer,M::integer,CB0::Vector(datatype=integer[4])) local i::integer,j::integer; CB0[1]:=1:

>	for i from 1 to N do CB0[i+1]:=CB0[i]+2:od:

>	for j from 2 to M-1 do i:=1:CB0[i+(j-1)N+1]:=CB0[i+(j-1)N]+2: for i from 2 to N-1 do CB0[i+(j-1)N+1]:=CB0[i+(j-1)N]+5:od: i:=N:CB0[i+(j-1)N+1]:=CB0[i+(j-1)N]+2:od:

>	for i from 1 to N do CB0[NM-N+1+i]:=CB0[NM-N+i]+2:od: end proc:

>	CBproc:=Compiler:-Compile(CBproc):

>

Rproc:=proc(N::integer,M::integer,R0::Vector(datatype=integer[4]))
local i::integer,j::integer,i1::integer,count::integer;
count:=0:
j:=1:
for i from 1 to N do
i1:=i+(j-1)*N:
#ff[i1]:=-y[i1]+y[i1+N]:
count:=count+1:
R0[count]:=i1:
count:=count+1:
R0[count]:=i1+N:
#j00[i1,i1]:=1.0:j00[i1,i1+N]:=-1.00:
od:

for j from 2 to M-1 do
i:=1:
i1:=i+(j-1)*N:
count:=count+1:
R0[count]:=i1:
count:=count+1:
R0[count]:=i1+1:
#j00[i1,i1]:=1.0:j00[i1,i1+1]:=-1.00:
for i from 2 to N-1 do
i1:=i+(j-1)*N:
count:=count+1:
R0[count]:=i1-N:
count:=count+1:
R0[count]:=i1-1:
count:=count+1:
R0[count]:=i1:
count:=count+1:
R0[count]:=i1+1:
count:=count+1:
R0[count]:=i1+N:
#j00[i1,i1]:=4*Y0[i,j]+Y0[i,j+1]+Y0[i,j-1]+Y0[i+1,j]+Y0[i-1,j]+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:
i:=N:
i1:=i+(j-1)*N:
count:=count+1:
R0[count]:=i1-1:
count:=count+1:
R0[count]:=i1:
i1:=i+(j-1)*N:
#j00[i1,i1]:=1.0:j00[i1,i1-1]:=-1.00:
od:

j:=M:
for i from 1 to N do
i1:=i+(j-1)*N:
count:=count+1:
R0[count]:=i1-N:
count:=count+1:
R0[count]:=i1:
#j00[i1,i1]:=1.0:j00[i1,i1-N]:=-1.00:
#ff[i1]:=-y[i1]+y[i1-N]+v0*h:
od:

end proc:

>	Rproc:=Compiler:-Compile(Rproc):

>

ENO2:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,Nx::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:

>

WENO3:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,Nx::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:=0.1:
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):
alpha:=1.0:beta:=1.0:
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:

>

UpW:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,Nx::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],e1::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:=0.1:
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:

>	ENO2:=Compiler:-Compile(ENO2):

>	WENO3:=Compiler:-Compile(WENO3):

>	UpW:=Compiler:-Compile(UpW):

>	printlevel:=1:

>

>	pFactor := define_external("hw_PardisoFactor", ':-MAPLE', ':-LIB' = "linalg"): pSolve := define_external("hw_PardisoSolve", ':-MAPLE', ':-LIB' = "linalg"): pCSF := define_external("hw_SpCompressedSparseForm", ':-MAPLE', ':-LIB' = "linalg"):

>

#NN:=20;

>	ss:=proc(NN,delPhi,nn,tf,beta,e1) local tt,V,V2,V1,TT,dt,vv0,ii,nt,Nt,N,h,w,Nx,MM,M,Ny,ntot,Ntot,Nsp,Y00,Phi0,F0,ff,db,CB,CB0,Ymid,R0,x0,handle,Phiave; N:=NN+2:h:=1.0/NN:w:=h/2:Nx:=NN:

>	MM:=NN:M:=MM+2:Ny:=MM:ntot:=N(M):Ntot:=ntot:Nsp:=5(N-2)^2+8*(N-1):

>	Y00:=Matrix(1..N,1..M,datatype=float[8]):F0:=Matrix(1..N,1..M,datatype=float[8]):Phi0:=Vector(1..Ntot,datatype=float[8]):ff:=copy(Phi0):db:=copy(ff):CB0:=Vector(N*M+1,datatype=integer[4]):Ymid:=Matrix(1..N,1..M,datatype=float[8]):

>	y0proc(NN,Y00,beta,e1):

>	R0:=Vector(Nsp,datatype=integer[4]):x0:=Vector(Nsp,datatype=float[8]):

>	CBproc(N,M,CB0):

>	CB0:=convert(CB0,Vector,datatype=integer[8]):

>	Rproc(N,M,R0):

>	R0:=convert(R0,Vector,datatype=integer[8]):

>	evalf(Eqs11(N,M,Y00,F0,Phi0,evalf(0.1),ff));

>	XX(N,M,Y00,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):#handle := pFactor(CB, R, X, 11): pSolve(ff,db,handle):

>	PhiAdd(Ntot,Phi0,db):

>	V2[0]:=(Phi0[ntot-2N+N/2]/2+Phi0[ntot-2N+N/2+1]/2)+h/20.1; V1[0]:=Phi0[2N]; V[0]:=Phi0[ntot-N]+0.5h0.1;

>	TT[0]:=0;tt:=0; 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/nn,tf-tt);

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

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

>	Nt:=round(tf/dt); #dt:=tf/Nt;

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

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

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

>

while tt<tf do
delPhi(Y00,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd(Y00,Ymid,Phi0,F0,dt,N,M,e1);
Eqs11(N,M,Ymid,F0,Phi0,0.1,ff);
XX(N,M,Ymid,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):
pSolve(ff,db,handle):
PhiAdd(Ntot,Phi0,db):
delPhi(Ymid,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd2(Y00,Ymid,Phi0,F0,dt,N,M,e1);
Eqs11(N,M,Ymid,F0,Phi0,evalf(0.1),ff);
XX(N,M,Ymid,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):#handle := pFactor(CB, R, X, 11):
pSolve(ff,db,handle):
PhiAdd(Ntot,Phi0,db):
delPhi(Ymid,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd3(Y00,Ymid,Phi0,F0,dt,N,M,e1);
Eqs11(N,M,Y00,F0,Phi0,evalf(0.1),ff);
XX(N,M,Y00,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):#handle := pFactor(CB, R, X, 11):
pSolve(ff,db,handle):
PhiAdd(Ntot,Phi0,db):
ii:=ii+1:#print(i);
V2[ii]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*0.1;#print(ii,V[ii]);
#V1[ii]:=Phi0[ntot]+h/2*0.1;
V1[ii]:=Phi0[2*N];
V[ii]:=Phi0[ntot-N]+0.5*h*0.1;
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:=h/vv0/nn;
dt:=min(h/vv0/nn,tf-tt);#print(tt,ii,dt);
gc();
end:

>

V2[ii];

>	end proc: #time[real]()-t11;time()-t12;NN;

>	ssapp:=proc(NN,delPhi,nn,tf,beta,e1) local tt,V,V2,V1,TT,dt,vv0,ii,nt,Nt,N,h,w,Nx,MM,M,Ny,ntot,Ntot,Nsp,Y00,Phi0,F0,ff,db,CB,CB0,Ymid,R0,x0,handle,Phiave, Phimid,Phitemp; N:=NN+2:h:=1.0/NN:w:=h/2:Nx:=NN:

>	MM:=NN:M:=MM+2:Ny:=MM:ntot:=N(M):Ntot:=ntot:Nsp:=5(N-2)^2+8*(N-1):

>

Y00:=Matrix(1..N,1..M,datatype=float[8]):F0:=Matrix(1..N,1..M,datatype=float[8]):Phi0:=Vector(1..Ntot,datatype=float[8]):ff:=copy(Phi0):db:=copy(ff):CB0:=Vector(N*M+1,datatype=integer[4]):Ymid:=Matrix(1..N,1..M,datatype=float[8]):Phitemp:=Vector(1..Ntot,datatype=float[8]):Phimid:=Vector(1..Ntot,datatype=float[8]):

>	y0proc(NN,Y00,beta,e1):

>	R0:=Vector(Nsp,datatype=integer[4]):x0:=Vector(Nsp,datatype=float[8]):

>	CBproc(N,M,CB0):

>	CB0:=convert(CB0,Vector,datatype=integer[8]):

>	Rproc(N,M,R0):

>	R0:=convert(R0,Vector,datatype=integer[8]):

>	evalf(Eqs11(N,M,Y00,F0,Phi0,evalf(0.1),ff));

>	XX(N,M,Y00,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):#handle := pFactor(CB, R, X, 11): pSolve(ff,db,handle):

>	PhiAdd(Ntot,Phi0,db):

>	V2[0]:=(Phi0[ntot-2N+N/2]/2+Phi0[ntot-2N+N/2+1]/2)+h/20.1; V1[0]:=Phi0[2N]; V[0]:=Phi0[ntot-N]+0.5h0.1;

>	TT[0]:=0;tt:=0; 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:=h/vv0/nn;

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

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

>	Nt:=round(tf/dt); #dt:=tf/Nt;

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

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

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

>

while tt<tf do
PhiSwitch(Ntot,Phi0,Phitemp):
delPhi(Y00,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd(Y00,Ymid,Phi0,F0,dt,N,M,e1);
Eqs11(N,M,Ymid,F0,Phi0,0.1,ff);
XX(N,M,Ymid,Phi0,0.1,x0);handle := pFactor(CB0, R0, x0, 11):
pSolve(ff,db,handle):
PhiAdd(Ntot,Phi0,db):
MidPred(Ntot,Phi0,Phitemp,Phimid):
delPhi(Ymid,Phi0,F0,evalf(dt),N,M,Nx,0.1):
EFAdd2(Y00,Ymid,Phi0,F0,dt,N,M,e1);
delPhi(Ymid,Phimid,F0,evalf(dt),N,M,Nx,0.1):
EFAdd3(Y00,Ymid,Phimid,F0,dt,N,M,e1);
ii:=ii+1:#print(i);
V2[ii]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*0.1;#print(ii,V[ii]);
#V1[ii]:=Phi0[ntot]+h/2*0.1;
V1[ii]:=Phi0[2*N];
V[ii]:=Phi0[ntot-N]+0.5*h*0.1;
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:=h/vv0/nn;
dt:=min(h/vv0,tf-tt);#print(tt,ii,dt);
gc();
end:
V2[ii];
end proc:
#time[real]()-t11;time()-t12;NN;

>

A procedure is written to find the value of the potential at t =tf using different upwind methods usinb both SSP RK3 and SSP RK3a. Nt1 can be set to 8 to get results till 1280.

>

Nt1:=5;

(4)

>	for i from 1 to Nt1 do ntot:=10*2^(i-1): t11:=time[real]():t12:=time():ss1[i]:=ss(ntot,WENO3,1,2.0,0.5,1e-9): if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end: print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

(5)

>	for i from 1 to Nt1 do ntot:=10*2^(i-1): t11:=time[real]():t12:=time():ss1[i]:=ss(ntot,ENO2,1,2.0,0.5,1e-9): if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end: print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

(6)

>	for i from 1 to Nt1 do ntot:=10*2^(i-1): t11:=time[real]():t12:=time():ss1[i]:=ss(ntot,UpW,1,2.0,0.5,1e-9): if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end: print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

(7)

>

Nt1:=5;

(8)

>	for i from 1 to Nt1 do ntot:=10*2^(i-1): t11:=time[real]():t12:=time():ss1[i]:=ssapp(ntot,WENO3,1,2.0,0.5,1e-9): if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end: print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

(9)

>

>	for i from 1 to Nt1 do ntot:=10*2^(i-1): t11:=time[real]():t12:=time():ss1[i]:=ssapp(ntot,ENO2,1,2.0,0.5,1e-9): if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end: print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

(10)

>

>	for i from 1 to Nt1 do ntot:=10*2^(i-1): t11:=time[real]():t12:=time():ss1[i]:=ssapp(ntot,UpW,1,2.0,0.5,1e-9): if i=1 then err:=1: else err:=ss1[i]-ss1[i-1]:end: print(time()-t12,time[real]()-t11,ntot,err,ss1[i]):od:

(11)

>

Download Pardisocodetoupload.mw

any takers?...

Answered: Dr. Venkat Subramanian 431

March 07 2022

0 0

Just checking to see if this can be done. We are in the process of writing a paper, and this can help in speeding up the code shared with the community.

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

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

Ask a Question

Create a Post

Venkat Subramanian

431 Reputation

13 Badges

MaplePrimes Activity

These are replies submitted by

Please try and let me know...

Tridiagonal is easy...

Read again...

not an answer but see this thread...

I believe we are done!...

Faster code almost done with UMFPACK so ...

fixed this...

Updated code...

thanks...

bump up...

of course Linux is ok to start with...

Windows please...

2 of the 3 suggestions...

Excellent...

any takers?...

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