conversion from marching formula to ode15s (PDE and MOL)

3 次查看（过去 30 天）
Moritz
Moritz 2013-7-7
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)));
  4 个评论
显示 2更早的评论
Moritz
Moritz 2013-7-9
Actually i am not sure wheter it is a stiff equation or not. Probably i should just try ode45
Moritz
Moritz 2013-7-9
编辑：Moritz 2013-7-9
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by