conversion from marching formula to ode15s (PDE and MOL)

2 views (last 30 days)
Hello,
i tried a toolbox for a convcection diffusion PDE (phenomenological sedimentation model) ( toolbox works fast but it is an external fortran routine) writing it with a marching formula is very slow. So i would like to solve it with ode15s but I probably do not understand how to vectorize the PDE. Using ode15s i get an error message.
So i do have two Questions and would appreciate any help:
A) Anyone got a similar error ? Any hint what i did wrong ?
B) what about a jacobian pattern ? i found an excample from mathworks where they provide it:
Thank you
Moritz
I do get follwing error message:
Error in ode15s (line 111) solver_name = 'ode15s';
Output argument "varargout{4}" (and maybe others) not assigned during call to "/usr/local/MATLAB/R2012b/toolbox/matlab/funfun/ode15s.m>ode15s".
Error in CallPhenSettlModel (line 72) [TOUT,YOUT,TE,YE,IE] = ode15s(f,TSPAN,v,OPTIONS,r0,rb,omega,u0,g,N,xcen,theta) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This post is organised as followed:
1) marching formula approach (upwind scheme) (working but slow in this way) 2) additional m files 3)try with ode15s
1) Marching Formula approach
function [v,FinalTime,xcen] = main
clear all, close all, clc
%%Initial Parameter
global u0 omega sigmaN k drho C g uinf r0 rb umax uc um
C=5; %Exponent in flux function
uc=0.07; %critical concentration, below this the problem will become hyperbolic
sigmaN=5.7; % material constant
k=9;% Exponent in flux function
umax=0.66; % material constant
u0=0.06; % Initial concentration
drho=1200; % density difference in kg/m^3
uinf=0.001; %material constant
r0=0.1; % radius at the liquor height
rb=0.6; % radius at the bottom of the tube
g=9.81; % earth acceleration m/s^2
omega=sqrt(100000/rb); % angular speed of centrifuge (assumption)
%%Assumed flux and diffusion function
uclose=linspace(0,1,100)';
au=afun(uclose); fu=ffun(uclose);
figure(1), subplot(3,1,1); plot(uclose,au,'r.'); ylabel('a');xlabel('uclose');
subplot(3,1,2); plot(uclose,fu,'r.');ylabel('fflux');xlabel('uclose');
%%The maximum of the flux function which is used in the numerical flux
[gm, im]=max(fu); um=uclose(im);
figure(1), subplot(3,1,2); hold on; plot(um,gm,'bo')
hold off
%%The integral A is calculated by the trapezoidal rule
% global Au
% Au=[0;cumsum(([au(1:length(uclose)-1);0]+[au(2:length(uclose));0])/2)];
Au=cumtrapz(uclose,afun(uclose));
figure(1),subplot(3,1,3); plot(uclose,Au,'r.');xlabel('uclose');ylabel('Au');
%%Main routine
% N = number of points to solve at
N = 50;
% range [a,b] to solve on
a = r0;
b = rb;
% compute spacing (here dx is constant)
dx=(b-a)/N;
% location of centerpoints
xcen = linspace(a+0.5*dx, b-0.5*dx, N);
% location of left cell interface
xl=xcen-0.5*dx;
%location of right interface
xr=xcen+0.5*dx;
% initial condition
u = u0.*ones(size(xcen))';
% u = linspace(0.01,1,size(xcen,2)).*ones(size(xcen,2),1)';
% compute time step
dt =0.9*dx/(b/a*(omega^2*uinf+2/dx*max(afun(uclose))));
% simulation time
StartTime = 0.;
FinalTime = 0.001;
% correct dt so that last time step ends at FinalTime
Nsteps = floor(FinalTime/dt);
% dt = FinalTime/Nsteps;
% tOut=StartTime:dt:FinalTime;
% Select Theta for limiter [0 2]
theta=0;
%%Main routine
timer=0;
% v=zeros(N,Nsteps);
v= u;%u(:,timer);
figure(2)
for n=1:Nsteps;
%pause(0.1)
timer=timer+1;
%calculate slope limiters at each center point including ghost cells
s=zeros(N,1);
for i=2:N-1
s(1)=0;
s(i)=minmod(theta.*(v(i)-v(i-1))./dx,(v(i+1)-v(i-1))./(2*dx),theta.*(v(i+1)-v(i))./dx);
s(N)=0;
end
%%Calculate extrapolated values at cell interfaces for interior
vL=v-dx/2*s;
vR=v+dx/2*s;
for xmesh=2:N-1
% using interpolation is faster
v(1)=v(1)- omega^2*dt/dx/g*(xl(1)*feo('ffun',v(1),v(2))) + ...
dt/dx^2*((interp1q(uclose,Au,v(2))-interp1q(uclose,Au,v(1))));
v(xmesh)=v(xmesh)- omega^2*dt/dx/g*((xmesh+0.5*dx)*feo('ffun',vR(xmesh),vL(xmesh+1))....
-(xmesh-0.5*dx)*feo('ffun',vR(xmesh-1),vL(xmesh))) ...
+ dt/dx^2*(interp1q(uclose,Au,v(xmesh+1))-interp1q(uclose,Au,v(xmesh))- ...
IntA(v(xmesh))-IntA(v(xmesh-1)));
v(N)=v(N1)+ omega^2*dt/dx/g*(xl(end)*feo('ffun',v(N-1),v(N))) - ...
dt/dx^2*(interp1q(uclose,Au,v(N))-interp1q(uclose,Au,v(N-1)));
% v=vn;
end
% This part is just for visual conrol during running
figure(2)
plot(xcen,v,'r.')
ylim([0 1])
drawnow
figure(3)
plot(xcen,s,'b.')
end figure plot(xcen,v,'r.') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2) additional m-files
function fEo = feo(g,u,v)
%The Engquist-Osher numerical flux function is implemented as
global um
fgu=feval(g,u); fgum=feval(g,um); fgv=feval(g,v); fguvum=fgu+fgv-fgum;
lu=u<=um; lv=v<=um;
fEo= lu.*lv.*fgu + lu.*~lv.*fguvum + ~lu.*lv.*fgum + ~lu.*~lv.*fgv;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = ffun(u)
global uinf umax C
f=uinf.*u.*(1-u/umax).^C.*(u>=0).*(u<=umax);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a= afun(u)
%diffusion coefficient
global sigmaN uc C k umax drho g uinf
a=sigmaN.*uinf./(drho*g).*(1-u).^C.*(u/uc).^(k-1).*(u>=uc).*(u<umax);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=IntA(u)
global uc umax
ucl=0:u/10:u;
A=trapz(ucl,afun(ucl)).*(u>=uc).*(u<umax).*(u>0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mm = minmod(a,b,c)
mm = zeros(size(a))';
mm = (a.*b.*c>0).*min([a,b,c]);
mm =mm + (a.*b.*c<0).*max([a,b,c]);
3) try to use ode15s
clear all, close all, clc
% This skript sets all parameters and calls ode15s
%%Initial Parameter
global u0 omega sigmaN k drho C g uinf r0 rb umax uc um
C=5; %Exponent in flux function
uc=0.07; %critical concentration, below this the problem will become hyperbolic
sigmaN=5.7; % material constant
k=9;% Exponent in flux function
umax=0.66; % material constant
u0=0.06; % Initial concentration
drho=1200; % density difference in kg/m^3
uinf=0.001; %material constant
r0=0.1; % radius at the liquor height
rb=0.6; % radius at the bottom of the tube
g=9.81; % earth acceleration m/s^2
omega=sqrt(100000/rb); % angular speed of centrifuge (assumption)
%%Assumed flux and diffusion function
uclose=linspace(0,1,100)';
au=afun(uclose); fu=ffun(uclose);
figure(1), subplot(3,1,1); plot(uclose,au,'r.'); ylabel('a');xlabel('uclose');
subplot(3,1,2); plot(uclose,fu,'r.');ylabel('fflux');xlabel('uclose');
%%The maximum of the flux function which is used in the numerical flux
[gm, im]=max(fu); um=uclose(im);
figure(1), subplot(3,1,2); hold on; plot(um,gm,'bo')
hold off
%%The integral A is calculated by the trapezoidal rule
% global Au
% Au=[0;cumsum(([au(1:length(uclose)-1);0]+[au(2:length(uclose));0])/2)];
Au=cumtrapz(uclose,afun(uclose));
figure(1),subplot(3,1,3); plot(uclose,Au,'r.');xlabel('uclose');ylabel('Au');
%%Main routine
% N = number of points to solve at
N = 50;
% range [a,b] to solve on
a = r0;
b = rb;
% compute spacing (here dx is constant)
dx=(b-a)/N;
% location of centerpoints
xcen = linspace(a+0.5*dx, b-0.5*dx, N);
% location of left cell interface
xl=xcen-0.5*dx;
%location of right interface
xr=xcen+0.5*dx;
initial condition
u = u0.*ones(size(xcen))';
% u = linspace(0.01,1,size(xcen,2)).*ones(size(xcen,2),1)';
% compute time step
dt =0.9*dx/(b/a*(omega^2*uinf+2/dx*max(afun(uclose))));
% simulation time
StartTime = 0.;
FinalTime = 0.001;
% correct dt so that last time step ends at FinalTime
Nsteps = floor(FinalTime/dt);
% dt = FinalTime/Nsteps;
% tOut=StartTime:dt:FinalTime;
% Select Theta for limiter [0 2]
theta=0;
%%Main routine
timer=0;
% v=zeros(N,Nsteps);
v= u;%u(:,timer);
TSPAN=[0 ;1];
OPTIONS = [];%odeset('Vectorized','on','Jpattern',jpattern(N));
f=@PhenSettlModel;
[TOUT,YOUT,TE,YE,IE] = ode15s(f,TSPAN,v,OPTIONS,r0,rb,omega,u0,g,N,xcen,theta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dvdt = PhenSettlModel(t,v,r0,rb,omega,u0,g,N,xcen,theta)
% Start and end point of mesh
a = r0;
b = rb;
% compute spacing (here dx is constant)
dx=(b-a)/N;
% start vector
v=ones(1,N)*u0;
%Preallocating
dydt=zeros(N,1);
%calculate slope limiters at each center point including ghost cells
s=zeros(N,1)';
for i=2:N-1
s(i)=minmod(theta.*(v(i)-v(i-1))./dx,(v(i+1)-v(i-1))./(2*dx),theta.*(v(i+1)-v(i))./dx);
end
%%Calculate extrapolated values at cell interfaces for interior points
vL=v-dx/2*s;
vR=v+dx/2*s;
%%evaluate dvdt at boundary
i=1;
dvdt(i,:)=v(i)- omega^2/dx/g*((xcen(1)-0.5*dx)*feo('ffun',v(i),v(i+1))) + ...
1/dx^2*((IntA(v(i+1))-IntA(v(i))));
% evaluate dvdt at interior grid points
i=2:N-1;
dvdt(i,:)=v(i)- omega^2/dx/g.*((xcen(i)+0.5).*dx.*feo('ffun',vR(i),vL(i+1))-(xcen(i)-0.5*dx)...
.*feo('ffun',vR(i-1),vL(i))) ...
+ 1/dx^2.*(IntA(v(i+1))-IntA(v(i))- ...
(IntA(v(i))-IntA(v(i-1))));
% evaluate dvdt at second boundary
i=N; dvdt(i,:)=v(i)+ omega^2/dx/g*((xcen(end)-0.5*dx)*feo('ffun',v(i-1),v(i)))-1/dx^2*(IntA(v(i))-IntA(v(i-1)));
  4 Comments
Moritz
Moritz on 9 Jul 2013
Actually i am not sure wheter it is a stiff equation or not. Probably i should just try ode45
Moritz
Moritz on 9 Jul 2013
Edited: Moritz on 9 Jul 2013
well, probably ode45 is the right joice. No error message and really fast. But the solution is not correct. I will let you know if it works.

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!