{This script describes solute transport through three hydrogeologic layers.
II want to use it to simulate vertical salt water seepage.When I simulate the pressure or omega it gives realistic results
for both but when I couple both variables and equations it gives bad results. I think there is something wrong with the definition of p.
Maybe you can give a look at it. Besides this mathematical problem it gives a BAD NEIGHBOR FAULT as I ask you when I
simulate the layers with (real) little thickness. Hopefullt you can help me with this. Thanks. greetz Marc Zwaanswijk}
TITLE '2D solute flow across Lichtaard, NW Friesland, Netherlands '
SELECT
stages=10
contours=30

VARIABLES     { system variables }
pressure
omega

DEFINITIONS   { parameter definitions }
k=1e-10
mu=1e-3
rho_0=1000
gy=10
g=vector(0,10)
D=9.5e-10
w = 0.0357*erf(-y/5)
rho=rho_0*exp(0.69*omega)
p=-rho*gy*y	
q=-(k/mu)*(grad(pressure)+rho*g)	
J=-rho*D*grad(omega)			{diffusive flux}
fac=10**((stage-10)/3)


INITIAL VALUES
omega=w	{ define initial distribution of concentration }
pressure=p	{ define initial distribution of pressure }

EQUATIONS     { PDE's, one for each variable } {voor stationaire situatie}
omega:div(rho*omega*q*fac)+div(J)=0		
pressure:div(rho*q*fac)=0

	
BOUNDARIES    { The domain definition }
  REGION 1 "total area"
    START(0,0)

	value(omega)=w
	value(pressure)=p

	line to (0,-32)

	value(omega) =0.0357
	natural(pressure)=0

	line to (2250,-32)

	value(omega)=w
	value(pressure)=p

	line to (2250,0)

	value(pressure)=0
	
	
              line to (1650,0)
 	



	line to (1650,-1.5)

	value(pressure)=0
	
			
	line to (1400,-1.5)



	line to (1400,0)

	value(pressure)=0
	
			
	line to (1150,0)



	line to (1150,-1.5)

	value(pressure)=0
	
			
	line to (950,-1.5)



	line to (950,0)

	value(pressure)=0
	
			
	line to (800,0)



	line to (800,-1.5)

	value(pressure)=0
	
			
	line to (600,-1.5)



	line to (600,0)

	value(pressure)=0
	
			
	line to finish

REGION 2 "deklaag"
 k=2.3e-14

  START(0,0)

	line to (0,-2.0)

	line to (2250,-2.0)

	line to (2250,0)

   	line to (1650,0)

	line to (1650,-1.5)

	line to (1400,-1.5)

	line to (1400,0)

	line to (1150,0)

	line to (1150,-1.5)

	line to (950,-1.5)

	line to (950,0)

	line to (800,0)

	line to (800,-1.5)

	line to (600,-1.5)

	line to (600,0)

	line to finish

 REGION 3 "fijnzandig wvp"
k=1.2e-13

  START(0,-2.0)

	line to (0,-9.0)

	line to (2250,-9.0)

	line to (2250,-2.0)

	line to finish

 REGION 4 "grofzandig wvp"
 k=1.2e-10

  START(0,-9.0)

	line to (0,-32)

	line to (2250,-32)

	line to (2250,-9.0)

	line to finish

REGION 5 "bodem afgraving 1"
 k=7e-14

  START(600,-1.5)

	line to (600,-2.0)

	line to (800,-2.0)

	line to (800,-1.5)

	line to finish

REGION 6 "bodem afgraving 2"
 k=7e-14

  START(950,-1.5)

	line to (950,-2.0)

	line to (1150,-2.0)

	line to (1150,-1.5)

	line to finish

REGION 7 "bodem afgraving 3"
 k=7e-14

  START(1400,-1.5)

	line to (1400,-2.0)

	line to (1650,-2.0)

	line to (1650,-1.5)

	line to finish

FEATURE START "afgraving 1" (600,-2.0) LINE TO (800,-2.0)
FEATURE START "afgraving 2" (950,-2.0) LINE TO (1150,-2.0)
FEATURE START "afgraving 3" (1400,-2.0) LINE TO (1650,-2.0)


MONITORS

contour(pressure) painted zoom
contour(w) painted zoom
contour(omega) painted zoom
contour(p) painted zoom
vector(q)  zoom
elevation(omega) from (0,0) to (0,-32)
elevation(w) from (0,0) to (0,-32)

PLOTS
contour(pressure) painted zoom
contour(w) painted zoom
contour(omega) painted zoom
contour(p) painted zoom
vector(q)  zoom
elevation(omega) from (0,0) to (0,-32)
elevation(w) from (0,0) to (0,-32)

SUMMARY
EXPORT FILE="Integrals"
REPORT LINE_INTEGRAL(normal(-(k/mu)*(grad(pressure)+rho*g)),"afgraving 1")
REPORT LINE_INTEGRAL(normal(-(k/mu)*(grad(pressure)+rho*g)),"afgraving 2")
REPORT LINE_INTEGRAL(normal(-(k/mu)*(grad(pressure)+rho*g)),"afgraving 3")



END