Nonlinear Problem - Darcy Flow

Try adding
parameters.parse()
to the beginning of your program and run with parameter
--petsc.info

see
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscInitialize.html
for more options.

Hi Felix,

Thanks for the quick reply! So I've added parameters.parse() to the beginning of the program and passed the --petsc.info parameter. I get:

[0] PetscInitialize(): PETSc successfully started: number of processors = 1
[0] PetscGetHostName(): Rejecting domainname, likely is NIS ubuntu-M505.(none)
[0] PetscInitialize(): Running on machine: ubuntu-M505
[0] SlepcInitialize(): SLEPc successfully started
[0] PetscCommDuplicate(): Duplicating a communicator 2140416 174653792 max tags = 2147483647
[0] PetscCommDuplicate(): returning tag 2147483647
[0] Petsc_DelComm(): Deleting PETSc communicator imbedded in a user MPI_Comm 2140416
[0] PetscCommDestroy(): Deleting PETSc MPI_Comm 174653792
[0] Petsc_DelCounter(): Deleting counter data in an MPI_Comm 174653792
[0] Petsc_DelComm(): Deleting PETSc communicator imbedded in a user MPI_Comm 174653792
[0] PetscCommDuplicate(): Duplicating a communicator 2140416 174653496 max tags = 2147483647
[0] PetscCommDuplicate(): returning tag 2147483647
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483646
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483645
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483644
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483643
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483642
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483641
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483640
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483639
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483638
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483637
t = 0
Solving nonlinear variational problem.
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483636
[0] PetscCommDuplicate(): returning tag 2147483635
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483634
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483633
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483632
  Applying boundary conditions to linear system.
  Applying boundary conditions to linear system.
  Applying boundary conditions to linear system.
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483631
  Assembling matrix over cells [=====> ] 15.1%
  Assembling matrix over cells [===========> ] 29.4%
  Assembling matrix over cells [===================> ] 51.0%
  Assembling matrix over cells [===========================> ] 72.6%
  Assembling matrix over cells [===================================> ] 94.1%
  Assembling matrix over cells [======================================] 100.0%
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2387 X 2387; storage space: 0 unneeded,41899 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 26
[0] Mat_CheckInode(): Found 2385 nodes out of 2387 rows. Not using Inode routines
  Applying boundary conditions to linear system.
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483630
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2387 X 2387; storage space: 0 unneeded,41899 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 26
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2387 X 2387; storage space: 0 unneeded,41899 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 26
  Applying boundary conditions to linear system.
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483629
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2387 X 2387; storage space: 0 unneeded,41899 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 26
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2387 X 2387; storage space: 0 unneeded,41899 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 26
  Applying boundary conditions to linear system.
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483628
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2387 X 2387; storage space: 0 unneeded,41899 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 26
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2387 X 2387; storage space: 0 unneeded,41899 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 26
  Solving linear system of size 2387 x 2387 (PETSc Krylov solver).
[0] PetscCommDuplicate(): returning tag 2147483627
  PETSc Krylov solver starting to solve 2387 x 2387 system.
[0] PetscCommDuplicate(): returning tag 2147483626
[0] PetscCommDuplicate(): returning tag 2147483625
[0] PetscCommDuplicate(): returning tag 2147483624
[0] PCSetUp(): Setting up new PC
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483623
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483622
[0] PetscCommDuplicate(): returning tag 2147483621
[0] PetscCommDuplicate(): Using internal PETSc communicator 2140416 174653496
[0] PetscCommDuplicate(): returning tag 2147483620

and it just hangs there. Do you think this has to do with memory issues or an error in my equations? Thanks again.

- Ben

Sorry, I really don't know.
Try adding monitor_convergence = true to the parameters of the Krylov
solver.
There is also a petsc option for debugging memory, but probably the cg
does not converge. The maximum_iterations is set to 10000 by default.
This could indicate an error in the equations.

Simplify it till you get a working system. Then add back complexity

  * Make it a stationary problem.
  * One place you divide with kappa, which is ~1e-13. That
    might create an ill-conditioned problem.
  * Set all parameters to 1 :)
  * Is the system symetric? Only then can you use CG. GMRES
    is a much more robust iterative solver.

Johan

On 11/30/2012 06:06 PM, Benjamin Mayersohn wrote:
> New question #215682 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/215682
>
> Dear all,
>
> I'm using FEniCS to model a PDE that couples flow through a porous medium with solute transport. Here are several PDF's that explain the equations/boundary conditions, show the geometry, display the results, etc.
>
> https://dl.dropbox.com/u/56831547/Darcy/2006WR004918.pdf
> https://dl.dropbox.com/u/56831547/Darcy/models.ssf.buoyancy_darcy_elder.pdf
>
> The second PDF describes how the problem is implemented in COMSOL, a commercial finite element program, and that's where I got the original idea. This one's the easier of the two to read through quickly.
>
> I've tried implementing the problem using the NonlinearVariationalSolver, as described in the "Nonlinear Problems" section in the FEniCS documentation. I'm initially setting dt=Constant(0.0) to make sure things are stable before advancing further in time. Here is my code:
>
> from dolfin import *
>
> # PARAMETERS (SI units)
> length = 300 # length of geometry
> height = 150 # height of geometry
> c_0 = 0.0 # initial concentration
> c_br = 1.0 # brine concentration
> theta = 0.1 # porosity
> rho_b = 1400.0 # bulk density
> rho_s = 1200.0 # brine density
> rho_0 = 1000.0 # pristine water density
> D_L = 3.56e-6 # molecular diffusion
> kappa = 4.9346e-13 # permeability
> beta = (rho_s-rho_0)/(c_br-c_0) # increase in density due to difference in salt concentration
> g_const = 9.81 # gravitational acceleration
> mu = .001 # dynamic viscosity
> p_0 = 0.0 # reference pressure
>
> # MESH
> scale = 10
> xnum = length/scale
> ynum = height/scale
>
> mesh = Rectangle(0,0,length,height,xnum,ynum,'left')
> n = FacetNormal(mesh)
>
> # FUNCTION SPACE
>
> # We will have two spaces, V1 and V2
> # P1: Solute concentration (discretized linearly)
> # P2: Pressure (discretized quadratically)
> P1 = FunctionSpace(mesh, "Lagrange", 1)
> P2 = FunctionSpace(mesh, "Lagrange", 2)
>
> # Combine into mixed function space ME
> ME = P1*P2
>
> # q, d are the test functions for p and c (respectively)
> du = TrialFunction(ME)
> d,q = TestFunctions(ME)
>
> u = Function(ME)
> u_prev = Function(ME) # Solution from previous time-step
>
> dc, dp = split(du)
> c, p = split(u)
> c_prev,p_prev = split(u_prev)
>
> # VARIABLES (PART 1)
> D = Expression('x[1]-height',height=height,cell=triangle) # elevation
>
> rho = rho_0 + beta*c # density
> rho_prev = rho_0 + beta*c_prev
> v = -mu/kappa*(grad(p)+rho*g_const*grad(D)) # Darcy velocity
> v_prev = -mu/kappa*(grad(p_prev)+rho_prev*g_const*grad(D))
>
>
> # initial pressure
> p_init_e = Expression('rho_0*g_const*(x[1]-height)',rho_0=rho_0,g_const=g_const,height=height)
>
> # BOUNDARY CONDITIONS
>
> # Dirchlet conditions
>
> def cbottom(x):
> return x[1] < DOLFIN_EPS
>
> cbrine = compile_subdomains("x[1]>height-DOLFIN_EPS && x[0]>=length/2. && x[0]<=length")
>
> # We have a reference pressure p_0 in the upper left hand corner
> # Enforce as a pointwise constraint
> def corner(x):
> return x[1]>height-DOLFIN_EPS and x[0]<DOLFIN_EPS
>
> bcs = ([DirichletBC(ME.sub(0),c_0,cbottom),
> DirichletBC(ME.sub(0),c_br,cbrine),
> DirichletBC(ME.sub(1),p_0,corner,method="pointwise")])
>
> # INITIAL CONDITIONS
> class InitialConditions(Expression):
> def __init__(self, c0, p0):
> self.c0 = c0
> self.p0 = p0
>
> def eval(self, value, x):
> value[0] = self.c0(x)
> value[1] = self.p0(x)
>
> def value_shape(self):
> return (2,)
>
> c_0 = Constant(c_0)
> c_init = interpolate(c_0,P1)
>
> p_init = interpolate(p_init_e,P2)
>
> u_init = InitialConditions(c_init,p_init)
> u.interpolate(u_init)
> u_prev.interpolate(u_init) # u_init no longer needed now
>
> # TIME-STEPPING
>
> # We'll use Crank-Nicholson
> # Let gamma be the evolution parameter
> gamma = 0.5
> p_mid = (1.0-gamma)*p_prev + gamma*p
> c_mid = (1.0-gamma)*c_prev + gamma*c
> rho_mid = (1.0-gamma)*rho_prev + gamma*rho
> v_mid = (1.0-gamma)*v_prev + gamma*v
>
> # WEAK FORMS
> dt = Constant(0.0) # Initially we plot the stationary case
>
> # For Darcy flow
> F0 = rho*(c-c_prev)*d*dx - dt/theta*(c_mid*div(rho_mid*d*v_mid))*dx+D_L*dt/2.0*inner(rho_prev*grad(c_prev)+rho*grad(c),grad(d))*dx
>
> # For solute transport
> F1 = q*(rho/beta+c)*div(rho*v)*dx + inner(rho*(c*v-theta*D_L*grad(c)),grad(q))*dx
>
> F = F0 + F1
> J = derivative(F,u,du)
>
> # SOLVE!
>
> # t = current time, T = end time
> t = 0
> T = 20*365*24*60*60
>
> while (t<=0):
> print "t =",t
>
> problem = NonlinearVariationalProblem(F, u, bcs, J)
> solver = NonlinearVariationalSolver(problem)
>
> prm = solver.parameters
> #info(prm, True)
> #print list_linear_solver_methods()
> #print list_krylov_solver_preconditioners()
>
> prm['linear_solver'] = 'cg'
> prm['preconditioner'] = 'ilu'
> prm['krylov_solver']['absolute_tolerance'] = 1E-3
> prm['krylov_solver']['relative_tolerance'] = 1E-6
> prm['krylov_solver']['maximum_iterations'] = 10
> prm['krylov_solver']['preconditioner']['ilu']['fill_level'] = 0
> set_log_level(PROGRESS)
>
> solver.solve()
> end()
>
> # Set dt to the Courant number
> # dt = h/u
> # dt = mesh.hmin()/u.max();
>
> t += dt
> u_prev.assign(u)
> #plot(p)
>
> But when I try solving the problem, I get:
>
> t = 0
> Solving nonlinear variational problem.
> Applying boundary conditions to linear system.
> Applying boundary conditions to linear system.
> Applying boundary conditions to linear system.
> Assembling matrix over cells [=====> ] 15.1%
> Assembling matrix over cells [======================> ] 58.0%
> Assembling matrix over cells [===========================> ] 72.3%
> Assembling matrix over cells [===================================> ] 93.9%
> Assembling matrix over cells [======================================] 100.0%
> Applying boundary conditions to linear system.
> Applying boundary conditions to linear system.
> Applying boundary conditions to linear system.
> Solving linear system of size 2387 x 2387 (PETSc Krylov solver).
> PETSc Krylov solver starting to solve 2387 x 2387 system.
>
> and from there it simply hangs. Do you have any suggestions as how I can determine what's wrong (i.e. get more output), and if there's anything else I'm missing. Thanks.
>
> - Ben
>