$onText
Stationary flow of an incompressible fluid in a rectangular area in the
presence of an obstacle.

Adapted from:
McKinney, D.C., Savitsky, A.G., Basic optimization models for water and
energy management. June 1999 (revision 6, February 2003).
http://www.ce.utexas.edu/prof/mckynney/ce385d/papers/GAMS-Tutorial.pdf
$offText

set X /u1*u20/;
set Y /u1*u20/;

* Determination of zone for water movement equation
Set vyside(X,Y);
Set vxside(X,Y);
    vxside(X,Y)                  = yes;
    vxside(X,Y)$(ord(X)=1)       = no;
    vxside(X,Y)$(ord(X)=card(X)) = no;
    vxside(X,Y)$(ord(Y)=1)       = no;
    vxside(X,Y)$(ord(Y)=card(Y)) = no;
    vxside('u10','u10')            = no;
    vxside('u10','u11')            = no;
    vyside(X,Y) = vxside(X,Y);

* Parameters
scalar dx  step space in x direction   /1/;
scalar dy  step space in y direction   /1/;
scalar r   density of the fluid        /1000/;

parameter m(X,Y) kinematic viscosity ;
          m(X,Y) := 0.5;

variables
  obj        objective variable
  D(X,Y)     error
  P(X,Y)     pressure
  Vx(X,Y)    x-direction velocity
  Vy(X,Y)    y-direction velocity  ;

* Variable bounds and initialization
  Vx.l(X,Y)     =  0.0;
  Vy.l(X,Y)     =  0.0;
  Vy.l(x,y)$(vxside(x,y)) = 0.0;

  D.lo(X,Y)     = 0.0;
  D.up(X,Y)     = 10.0;

  P.up(X,Y)     =  20000;
  P.lo(X,Y)     = -20000;
  P.l(X,Y)      =  0.0;

*Boundary conditions
  Vx.fx('u1',Y)   = 0.5;
  Vx.fx('u20',Y)  = 0.5;
  Vx.fx(X,'u1')   = 0;
  Vx.fx(X,'u20')  = 0;
  Vy.fx('u1',Y)   = 0;
  Vy.fx('u20',Y)  = 0;
  Vy.fx(X,'u1')   = 0;
  Vy.fx(X,'u20')  = 0;

* Obstacle description
  vx.fx('u10','u10') = 0;
  vx.fx('u10','u11') = 0;

Equations
     For_Vx(X,Y)
     For_Vy(X,Y)
     Div_Vxy(X,Y)
     eobj  objective function ;

For_Vx(X,Y)$(vxside(X,Y))..
   (P(X+1,Y)-P(X,Y))/(r*dx) =e=
     m(x,Y)*((Vx(X+1,Y)-2*Vx(X,Y)+Vx(x-1,Y))/(dx*dy)
    +        (Vx(X,Y+1)-2*Vx(X,Y)+Vx(X,Y-1))/(dx*dy));

For_Vy(X,Y)$(vxside(X,Y))..
   (P(X,Y+1)-P(X,Y))/(r*dy) =e=
     m(X,Y)*((Vy(X+1,Y)-2*Vy(X,Y)+Vy(X-1,Y))/(dx*dy)
    +        (Vy(X,Y+1)-2*Vy(X,Y)+Vy(X,Y-1))/(dy*dx));
*--
*For_Vx(X,Y)$(Vxside(X,Y))..
*    Vx(X,Y)*(Vx(X+1,Y)-Vx(X-1,Y))/(2*dx) +
*    0.25*(Vy(X+1,Y-1)+Vy(X+1,Y)+Vy(X,Y-1)+Vy(X,Y)) *
*         (Vx(X,Y+1)-Vx(X,Y-1))/(2*dy) +
*         (P(X+1,Y)-P(X,Y))/(r*dx)
*    =e=
*    m(X,Y)*((Vx(X+1,Y)-2*Vx(X,Y)+Vx(X-1,Y))/(dx*dx) +
*            (Vx(X,Y+1)-2*Vx(X,Y)+Vx(X,Y-1))/(dy*dy));
*
*For_Vy(X,Y)$(Vyside(X,Y))..
*    0.25*(Vx(X-1,Y+1)+Vx(X-1,Y)+Vx(X,Y+1)+Vx(X,Y)) *
*         (Vy(X+1,Y)-Vy(X-1,Y))/(2*dy) +
*         Vy(X,Y)*(Vy(X,Y+1)-Vy(X,Y-1))/(2*dy) +
*         (P(X,Y+1)-P(X,Y))/(r*dy)
*    =e=
*    m(X,Y)*((Vy(X+1,Y)-2*Vy(X,Y)+Vy(X-1,Y))/(dx*dx) +
*            (Vy(X,Y+1)-2*Vy(X,Y)+Vy(X,Y-1))/(dy*dy));
*--
Div_Vxy(X,Y)$((ord(X) > 1) $ (ord(Y) > 1))..
    (Vx(X,Y)-Vx(X-1,Y))/dx + (Vy(X,Y)-Vy(X,Y-1))/dy =e=
         D(X,Y);

eobj.. obj =e= SUM((X,Y),(D(X,Y)*D(X,Y)));

Model flowobs /all/;

Solve flowobs using nlp minimizing obj;

$ifThenI x%mode%==xbook
* Put the solution
  file res1 /pressure.dat/
  put res1;
  loop(X, loop(Y, put P.l(x,y):10:4; ); put /;); put /;

  file res2 / vx.dat/
  put res2;
  loop(X, loop(Y, put Vx.l(x,y):9:5; ); put /;); put /;

  file res3 / vy.dat/
  put res3;
  loop(X, loop(Y, put Vy.l(x,y):5:1; ); put /;); put /;

  file res4 /err.dat/
  put res4
  loop(X, loop(Y, put D.l(x,y):5:1; ); put /;); put /;
$endIf

* end flowobs
*************************************************** November 19, 2008