sundialsTB > cvodes > cvbx.m
cvbx

PURPOSE

CVBX - CVODES example problem (serial, band)

SYNOPSIS

This is a script file.

DESCRIPTION

CVBX - CVODES example problem (serial, band)
   The following is a simple example problem with a banded Jacobian,
   with the program for its solution by CVODES.
   The problem is the semi-discrete form of the advection-diffusion
   equation in 2-D:
      du/dt = d^2 u / dx^2 + .5 du/dx + d^2 u / dy^2
   on the rectangle 0 <= x <= 2, 0 <= y <= 1, and the time
   interval 0 <= t <= 1. Homogeneous Dirichlet boundary conditions
   are posed, and the initial condition is
      u(x,y,t=0) = x(2-x)y(1-y)exp(5xy).
   The PDE is discretized on a uniform MX+2 by MY+2 grid with
   central differencing, and with boundary values eliminated,
   leaving an ODE system of size NEQ = MX*MY.
   This program solves the problem with the BDF method, Newton
   iteration with the CVBAND band linear solver, and a user-supplied
   Jacobian routine.
   It uses scalar relative and absolute tolerances.
   Output is printed at t = .1, .2, ..., 1.
   Run statistics (optional outputs) are printed at the end.

   See also: cvbx_f, cvbx_q, cvbx_J

CROSS-REFERENCE INFORMATION

This function calls:
  • CVode CVode integrates the ODE.
  • CVodeFree CVodeFree deallocates memory for the CVODES solver.
  • CVodeGetStats CVodeGetStats returns run statistics for the CVODES solver.
  • CVodeMalloc CVodeMalloc allocates and initializes memory for CVODES.
  • CVodeMonitor CVodeMonitor is the default CVODES monitoring function.
  • CVodeSetOptions CVodeSetOptions creates an options structure for CVODES.
  • cvbx_J CVBX_J - Jacobian functino for the CVBX example problem.
  • cvbx_f CVBX_F - RHS function for the CVBX example problem
  • cvbx_q CVBX_Q - quadrature function for the CVBX example problem.
This function is called by:

SOURCE CODE

0001 %CVBX - CVODES example problem (serial, band)
0002 %   The following is a simple example problem with a banded Jacobian,
0003 %   with the program for its solution by CVODES.
0004 %   The problem is the semi-discrete form of the advection-diffusion
0005 %   equation in 2-D:
0006 %      du/dt = d^2 u / dx^2 + .5 du/dx + d^2 u / dy^2
0007 %   on the rectangle 0 <= x <= 2, 0 <= y <= 1, and the time
0008 %   interval 0 <= t <= 1. Homogeneous Dirichlet boundary conditions
0009 %   are posed, and the initial condition is
0010 %      u(x,y,t=0) = x(2-x)y(1-y)exp(5xy).
0011 %   The PDE is discretized on a uniform MX+2 by MY+2 grid with
0012 %   central differencing, and with boundary values eliminated,
0013 %   leaving an ODE system of size NEQ = MX*MY.
0014 %   This program solves the problem with the BDF method, Newton
0015 %   iteration with the CVBAND band linear solver, and a user-supplied
0016 %   Jacobian routine.
0017 %   It uses scalar relative and absolute tolerances.
0018 %   Output is printed at t = .1, .2, ..., 1.
0019 %   Run statistics (optional outputs) are printed at the end.
0020 %
0021 %   See also: cvbx_f, cvbx_q, cvbx_J
0022 
0023 % Radu Serban <radu@llnl.gov>
0024 % Copyright (c) 2005, The Regents of the University of California.
0025 % $Revision: 1.3 $Date: 2006/03/07 01:19:54 $
0026 
0027 xmax = 2.0;
0028 ymax = 1.0;
0029 mx = 10;
0030 my = 5;
0031 
0032 rtol = 0.0;
0033 atol = 1.0e-5;
0034 
0035 t0 = 0.0;
0036 dtout = 0.1;
0037 nout = 10;
0038 
0039 dx = xmax/(mx+1);
0040 dy =  ymax/(my+1);
0041 
0042 % Problem data structure
0043 data.xmax = xmax;
0044 data.ymax = ymax;
0045 data.mx = mx;
0046 data.my = my;
0047 data.dx = dx;
0048 data.dy = dy;
0049 data.hdcoef = 1.0/dx^2;
0050 data.hacoef = 0.5/(2.0*dx);
0051 data.vdcoef = 1.0/dy^2;
0052 
0053 % initial conditions for states
0054 t = t0;
0055 u = zeros(mx*my,1);
0056 for j = 1:my
0057   y = j * data.dy;
0058   for i = 1:mx
0059     x = i * data.dx;
0060     u(j + (i-1)*my) = x*(xmax-x)*y*(ymax-y)*exp(5.0*x*y);
0061   end
0062 end
0063 
0064 % Initial condition for quadrature variable
0065 q = 0.0;
0066 
0067 options = CVodeSetOptions('RelTol',rtol, 'AbsTol',atol);
0068 options = CVodeSetOptions(options,...
0069                           'LinearSolver','Band',...
0070                           'JacobianFn',@cvbx_J,...
0071                           'UpperBwidth',my,...
0072                           'LowerBwidth',my);
0073 options = CVodeSetOptions(options,...
0074                           'Quadratures', 'on',...
0075                           'QuadRhsFn',@cvbx_q,...
0076                           'QuadInitCond', q);
0077 
0078 
0079 
0080 mondata.grph = false;
0081 options = CVodeSetOptions(options,...
0082                           'MonitorFn',@CVodeMonitor,...
0083                           'MonitorData',mondata);
0084 
0085 CVodeMalloc(@cvbx_f, t, u, options, data);
0086 
0087 
0088 ff=figure;
0089 hold on;
0090 box
0091 
0092 umax = norm(u,'inf');
0093 uavg = cvbx_q(t,u,data);
0094 fprintf('At t = %f   max.norm(u) = %e\n',t, umax);
0095 
0096 for i = 1:nout
0097 
0098   t_old = t;
0099   uavg_old = uavg;
0100 
0101   tout = t + dtout;
0102   [status,t,u, q] = CVode(tout, 'Normal');
0103 
0104   if status ~= 0
0105     return
0106   end
0107   
0108   uavg = cvbx_q(t,u,data);
0109   umax = norm(u,'inf');
0110   fprintf('At t = %f   max.norm(u) = %e\n',tout, umax);
0111 
0112   figure(ff);
0113   plot([t_old t],[uavg_old uavg]);
0114   plot([t0 tout], [q q]/(tout-t0), 'r-');
0115   plot([tout tout], [0 q]/(tout-t0), 'r-');
0116   
0117 end
0118   
0119 si= CVodeGetStats
0120 
0121 CVodeFree;
Generated on Thu 05-Oct-2006 15:56:58 by m2html © 2003