function pclaw_ex_obstacle_disode
%PCLAW_EX_OBSTACLE_DISODE
%   Example for PARTICLECLAW, a numerical scheme for scalar
%   conservation and balance laws
%       u_t+(f(u))_x = g(x,u)
%   in one space dimension.
%   
%   Example: Flow over obstacle
%       f(u) = u^2/2
%       g(x,u) = phi(x)*u
%   with Riemann initial data,
%   where phi(x) = pi*sin(pi*x)*(4.5<=x&x<=5.5)
%   is the gradient of the ground profile.
%
%   The source g is discontinuous. Solving the characteristic
%   equations as black box, as done in PCLAW_EX_OBSTACLE,
%   would yield only first order accuracy. In contrast, here
%   a problem specific ODE solver is provided.
%
%   Copyright (c) 2008 Benjamin Seibold and Yossi Farjoun
%   http://math.mit.edu/~seibold/research/particleclaw
%   http://arxiv.org/abs/0809.0726

%===============================================================================
parameters = struct(...
'name','Flow over obstacle',...  % name of example
'f',@f,'f1',@f1,'f2',@f2,...     % flux function and derivatives (below)
'ic',@ic,...                     % initial conditions (below)
'xbox',[0 12],...                % computational domain
'ubox',[-.1 3.1],...             % axis for function values (for plotting only)
'd',[0 1.25e-1 inf],...          % [d_min d_init d_max]
'tfinal',6,...                   % final time
'dtmax',1e-3,...                 % maximum time step (only relevant for sources)
'steps',250,...                  % number of output steps
'odestep',@odestep,...           % special method for characteristic ODE (below)
'flag_plot',1,...                % if yes, animation gets plotted in figure
'flag_sharpen',0,...             % if yes, shocks get sharpened upon output
'flag_save',0      );            % if yes, solution gets stored in data files
%===============================================================================
particleclaw(parameters)
%===============================================================================
function y = f(u),  y = u.^2/2;
function y = f1(u), y = u;
function y = f2(u), y = u*0+1;
function y = ic(x), y = 1+(x<1);

function [xn,un] = odestep(par,x,u,dt)
% Time step with implicit midpoint rule for the ODE x'=u,u'=phi(x)u.
% All points crossing a discontinuity in phi(x) are integrated
% in two steps: First up to the discontinuity, then boyond.
% input:  x,u    positions and values of points
%         dt     time step
% output: xn,un  positions and values after time step dt
xn = x; un = u;
if length(dt)==1, dt = x*0+dt; end
for ii = 1:5
    phin = phi((x+xn)/2);
    phi1n = phi1((x+xn)/2);
    % function F
    f1 = xn-dt/2.*(un+u)-x;
    f2 = un-dt.*phin.*(u+un)/2-u;
    % Jacobi matrix DF
    df11 = xn*0+1; df12 = -dt/2;
    df21 = -dt/2.*phi1n.*(u+un)/2; df22 = 1-dt/2.*phin;
    d = (df11.*df22-df12.*df21);
    xn = xn-(df22.*f1-df12.*f2)./d;
    un = un-(-df21.*f1+df11.*f2)./d;
end
% smooth integration over discontinuities
for xbar = [4.5 5.5];
    ind = find(x<xbar-1e-14&xn>xbar+1e-14);
    if any(ind) % find points that crossed a discontinuity
        phin = phi((x(ind)+xbar)/2);
        % go up to discontinuity
        xn(ind) = xbar;
        un(ind) = u(ind)+phin.*(xbar-x(ind));
        % do rest of step
        dt_loc = 1./(phin/2+u(ind)./(xbar-x(ind)));
        [xn(ind),un(ind)] = odestep(par,xn(ind),un(ind),dt(ind)-dt_loc);
    end
end

function y = phi(x)
% Derivative of ground profile.
y = pi*sin(pi*x).*(4.5<=x&x<=5.5);

function y = phi1(x)
% Second derivative of ground profile.
y = pi^2*cos(pi*x).*(4.5<=x&x<=5.5);

%===============================================================================
% Copyright (c) 2008 Benjamin Seibold and Yossi Farjoun
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction for non-commercial purposes, including
% without limitation the rights to use, copy, modify, merge, publish, and/or
% distribute copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software, and credit has to be
% given to the authors in publications that are in any form based on this
% Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
%===============================================================================