higher numerical precision

Please take a look at the code below. For drift-diffusion type problems my solution u is sometimes constant in certain regions of the mesh resulting in a zero gradient. This is problematic because I need to calculate

    J = n*u.dx(0)

where n has very big values and u.dx(0) has very small values. This results in J=0 which is unphysical, instead of the more accurate small number I want to calculate. I tried rescaling the problem and such but that did not help. Please see the code below that illustrates the problem. If Doping is too high (1e5) spikes start to appear in the current density Jn. Reduce Doping to 1e3 and the spikes disappear.

I am looking for suggestions on how to implement float128 precision or use a higher than double precision linear algebra library with dolfin. I know PETSc can be compiled with float128. Is it possible to use Epetra with float128? Any suggestion or experience with this type of problems is welcome.

    "Author: Chaffra Affouda, copyright 2014"
    from dolfin import *
    from dolfin.cpp.la import NewtonSolver

    parameters['linear_algebra_backend'] = 'PETSc'


    class DDProblem(NonlinearProblem):
      def __init__(self, F, J, bcs=None, form_compiler_parameters=None, ident_zeros=False):
        NonlinearProblem.__init__(self)
        self._F = F
        self._J = J
        self.bcs = bcs
        self.form_compiler_parameters = form_compiler_parameters
        self.ident_zeros = ident_zeros

      def F(self, b, x):
        assemble(self._F, tensor=b,
                 form_compiler_parameters=self.form_compiler_parameters)

        #Apply boundary conditions
        for bc in self.bcs:
            bc.apply(b, x)

      def J(self, A, x):


        assemble(self._J, tensor=A,
                 form_compiler_parameters=self.form_compiler_parameters)

        #Apply boundary conditions
        for bc in self.bcs:
            bc.apply(A)

        if self.ident_zeros:
            A.ident_zeros()


    class DDSolver(NewtonSolver):

      def __init__(self, du_max=1.1):
        self.du_max = du_max
        NewtonSolver.__init__(self)

      def update_solution(self, x, dx):
        self.damp_potentials(x, dx)
        NewtonSolver.update_solution(self, x, dx)

      def damp_potentials(self, x, dx):
        """This will strongly influence the convergence"""

        du = dx.array()
        dumax = self.du_max

        i = (du > dumax)
        du[i] = dumax#*(1 + np.log(du[i])/dumax)
        i = (du < -dumax)
        du[i] = -dumax#*(1 + np.log(-du[i])/dumax)
        dx[:] = du

    t1 = 1
    t2 = 1
    t3 = 0.1

    x0 = 0.0
    x1 = x0 + t1
    x2 = x1+t2
    x3 = x2+t3

    mesh = IntervalMesh(100,x0,x2)

    hmax = mesh.hmax()
    def inside_region1(x, on_boundary):
       return between(x[0],(x0,x1+hmax))
    def inside_region2(x, on_boundary):
      return between(x[0],(x1,x2+hmax))

    cell_domains = CellFunction('size_t',mesh)
    region1 = AutoSubDomain(inside_region1)
    region2 = AutoSubDomain(inside_region2)

    region1.mark(cell_domains,1)
    region2.mark(cell_domains,2)

    Q = FunctionSpace(mesh,'CG',1)
    DG1 = FunctionSpace(mesh,'DG',1)
    W = MixedFunctionSpace([Q,Q,Q])

    class Doping(Expression):
      def eval(self, v, x):
        if inside_region1(x,False):
            v[0] = -1e5
        elif inside_region2(x,False):
            v[0] = 1e5

    class Ni(Expression):
      def eval(self, v, x):
        if inside_region1(x,False):
            v[0] = 1e-5
        elif inside_region2(x,False):
            v[0] = 1e-5

    u_mixed = Function(W)
    v_mixed = TestFunction(W)
    du_mixed = TrialFunction(W)            
    u, un, up = split(u_mixed)
    v, vn, vp = split(v_mixed)
    du, dn, dp = split(du_mixed)

    epsilon = Constant(1.)
    T = Constant(1.)
    ni = Ni() #Constant(1.0) #Ni()
    d = Doping()
    n = ni*exp((u-un)/T)
    p = ni*exp((up-u)/T)
    q = 1.0
    rho = (p+d-n)*q
    mun = Constant(100.0)
    mup = Constant(100.0)

    tau_n = Constant(1.0)
    tau_p = Constant(1.0)
    R = (n*p - ni**2) + (n*p-ni**2)/(tau_n*(p+ni**2)+tau_p*(n+ni**2))

    dx = Measure("dx")[cell_domains]
    dS = Measure("dS")

    L1 = inner(epsilon*grad(u),grad(v))*dx(1) - rho*v*dx(1)
    L1 += inner(epsilon*grad(u),grad(v))*dx(2) - rho*v*dx(2)

    L2 = inner(mun*n*grad(un), grad(vn))*dx(1) - R*vn*dx(1)
    L2 += inner(mun*n*grad(un), grad(vn))*dx(2) - R*vn*dx(2)

    L3 = inner(mup*p*grad(up), grad(vp))*dx(1) + R*vp*dx(1)
    L3 += inner(mup*p*grad(up), grad(vp))*dx(2) + R*vp*dx(2)


    F = L1+L2+L3

    def inside_top_contact(x, on_boundary):
      return near(x[0],x2,)
    def inside_bottom_contact(x, on_boundary):
      return near(x[0],x0,)

    top_contact = AutoSubDomain(inside_top_contact)
    bottom_contact = AutoSubDomain(inside_bottom_contact)

    bias = Constant(1.0)

    J = derivative(F,u_mixed,du_mixed)
    F0 = L1
    J0 = derivative(F0,u,du)

    F = Form(F)
    J = Form(J)

    solver = DDSolver()
    solver.parameters['maximum_iterations'] = 100 
    solver.parameters['absolute_tolerance'] = 1e-10
    solver.parameters['relative_tolerance'] = 1e-30
    solver.parameters['relaxation_parameter'] = 1.0
    solver.parameters['error_on_nonconvergence'] = False
    solver.parameters['convergence_criterion'] = 'incremental'
    solver.parameters['linear_solver'] = 'lu'

    u_mixed.vector()[:] = 0.0
    problem_phi = DDProblem(F0, J0, bcs=[], ident_zeros=True)
    solver.solve(problem_phi, u_mixed.vector())
    plot(u, title='equilibrium potential',
         interactive=False)
    plt = plot(ln(abs(n))/ln(10), title="equilibrium electron density", interactive=False,          lc='1,0,0')

    #bc_bottom_phi = DirichletBC(W.sub(0),u0,bottom_contact)
    bc_bottom_n = DirichletBC(W.sub(1),0.0,bottom_contact)
    bc_bottom_p = DirichletBC(W.sub(2),0.0,bottom_contact)

    #bc_top_phi = DirichletBC(W.sub(0),u0+bias,top_contact)
    bc_top_n = DirichletBC(W.sub(1),bias,top_contact)
    bc_top_p = DirichletBC(W.sub(2),bias,top_contact)

    bcs_phi = []
    bcs_n  = [bc_bottom_n, bc_top_n]
    bcs_p = [bc_bottom_p, bc_top_p]

    problem_dd = DDProblem(F,J,bcs=bcs_phi+bcs_n+bcs_p,  ident_zeros=False)

    solver.parameters['relaxation_parameter'] = 1.0
    solver.du_max = 1.0
    solver.solve(problem_dd, u_mixed.vector())

    Jn = (-q*mun*n*un.dx(0))
    #Jp = (-q*mup*p*up.dx(0))
    #f = project(Jn+Jp,DG1)

    #plot(f,title='total current')

    plt= plot(ln(abs(Jn)+1e-10)/ln(10),title='Jn', lc='1,0,0')
    plt = plot(un, title='un',lc='1,0,0')

Comments (4)