conversion from marching formula to ode15s (PDE and MOL)

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

For more information about the PDE: http://www.ing-mat.udec.cl/~rburger/papers/bcijmp01.pdf*

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

[Marc]

Unfortunately, the code you posted above is very hard to follow due to the font changes. Maybe you can try and clean it up a bit and we can see if we can help.

You also mentioned a toolbox but failed to provide a reference for that. Can this be made available?

Since you are using or implying the use of method of lines, does your resulting PDEs become odes or a DAE?

[Moritz]

Thank you for the comment. I am working on it.

The toolbox is from NAG and i use the routine d03pf for it. I would like to write it in Matlab (or eventually in c++) in order to be compatible with other people.

[Moritz]

Actually i am not sure wheter it is a stiff equation or not. Probably i should just try ode45

[Moritz]

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.