So I am implementing PETSc/TAO's semi smooth solver to ensure non-negative solutions for anisotropic diffusion. Attached is the code.

Basically, the call-back functions require a Jacobian and a gradient/constraint evaluation routine. It begins around line 160. Run the code as:

python SS_ex2.py 30 <guess>

where <guess> is either use zeros as initial guess (0) or use the standard formulation solution (i.e., no bounds) as the initial guess (1). The implementation somewhat works if <guess> is 1 or an initial non-zero solution is provided. If no guess is provided, I get a funky solution where everywhere is zero except for the non-homogeneous DirichletBCs.

Also, even if I have an initial guess (1), the values of the BC aren't exactly accurate either (I should have a max value of 1 at the very center of the top face, but I get 0.983941)

Something tells me I am not implementing the BCs correctly. Can anyone tell me what I am possibly doing wrong here?

Thanks!

Justin