Python Function.assign方法代码示例

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
 def __init__(self, V):
     u = Function(V)
     u.assign(Constant('1.0'))
     self.one = u.vector()
     self.Mdiag = self.one.copy()
     self.Mdiag2 = self.one.copy()
     self.invMdiag = self.one.copy()

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
def _compute_time_errors(problem, method, mesh_sizes, Dt, plot_error=False):
    mesh_generator, solution, ProblemClass, cell_type = problem()
    # Translate data into FEniCS expressions.
    fenics_sol = Expression(smp.printing.ccode(solution['value']),
                            degree=solution['degree'],
                            t=0.0,
                            cell=cell_type
                            )
    # Compute the problem
    errors = {'theta': numpy.empty((len(mesh_sizes), len(Dt)))}
    # Create initial state.
    # Deepcopy the expression into theta0. Specify the cell to allow for
    # more involved operations with it (e.g., grad()).
    theta0 = Expression(fenics_sol.cppcode,
                        degree=solution['degree'],
                        t=0.0,
                        cell=cell_type
                        )
    for k, mesh_size in enumerate(mesh_sizes):
        mesh = mesh_generator(mesh_size)
        V = FunctionSpace(mesh, 'CG', 1)
        theta_approx = Function(V)
        theta0p = project(theta0, V)
        stepper = method(ProblemClass(V))
        if plot_error:
            error = Function(V)
        for j, dt in enumerate(Dt):
            # TODO We are facing a little bit of a problem here, being the
            # fact that the time stepper only accept elements from V as u0.
            # In principle, though, this isn't necessary or required. We
            # could allow for arbitrary expressions here, but then the API
            # would need changing for problem.lhs(t, u).
            # Think about this.
            stepper.step(theta_approx, theta0p,
                         0.0, dt,
                         tol=1.0e-12,
                         verbose=False
                         )
            fenics_sol.t = dt
            #
            # NOTE
            # When using errornorm(), it is quite likely to see a good part
            # of the error being due to the spatial discretization.  Some
            # analyses "get rid" of this effect by (sometimes implicitly)
            # projecting the exact solution onto the discrete function
            # space.
            errors['theta'][k][j] = errornorm(fenics_sol, theta_approx)
            if plot_error:
                error.assign(project(fenics_sol - theta_approx, V))
                plot(error, title='error (dt=%e)' % dt)
                interactive()
    return errors, stepper.name, stepper.order

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
 def assemble_solution(self, t):  # returns Womersley sol for time t
     if self.tc is not None:
         self.tc.start('assembleSol')
     sol = Function(self.solutionSpace)
     dofs2 = self.solutionSpace.sub(2).dofmap().dofs()  # gives field of indices corresponding to z axis
     sol.assign(Constant(("0.0", "0.0", "0.0")))  # QQ not needed
     sol.vector()[dofs2] += self.factor * self.bessel_parabolic.vector().array()  # parabolic part of sol
     for idx in range(8):  # add modes of Womersley sol
         sol.vector()[dofs2] += self.factor * cos(self.coefs_exp[idx] * pi * t) * self.bessel_real[idx].vector().array()
         sol.vector()[dofs2] += self.factor * -sin(self.coefs_exp[idx] * pi * t) * self.bessel_complex[idx].vector().array()
     if self.tc is not None:
         self.tc.end('assembleSol')
     return sol

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
def compute_velocity_correction(
    ui, p0, p1, u_bcs, rho, mu, dt, rotational_form, my_dx, tol, verbose
):
    """Compute the velocity correction according to

    .. math::

        U = u_0 - \\frac{dt}{\\rho} \\nabla (p_1-p_0).
    """
    W = ui.function_space()
    P = p1.function_space()

    u = TrialFunction(W)
    v = TestFunction(W)
    a3 = dot(u, v) * my_dx
    phi = Function(P)
    phi.assign(p1)
    if p0:
        phi -= p0
    if rotational_form:
        r = SpatialCoordinate(W.mesh())[0]
        div_ui = 1 / r * (r * ui[0]).dx(0) + ui[1].dx(1)
        phi += mu * div_ui
    L3 = dot(ui, v) * my_dx - dt / rho * (phi.dx(0) * v[0] + phi.dx(1) * v[1]) * my_dx
    u1 = Function(W)
    solve(
        a3 == L3,
        u1,
        bcs=u_bcs,
        solver_parameters={
            "linear_solver": "iterative",
            "symmetric": True,
            "preconditioner": "hypre_amg",
            "krylov_solver": {
                "relative_tolerance": tol,
                "absolute_tolerance": 0.0,
                "maximum_iterations": 100,
                "monitor_convergence": verbose,
            },
        },
    )
    # u = project(ui - k/rho * grad(phi), V)
    # div_u = 1/r * div(r*u)
    r = SpatialCoordinate(W.mesh())[0]
    div_u1 = 1.0 / r * (r * u1[0]).dx(0) + u1[1].dx(1)
    info("||u||_div = {:e}".format(sqrt(assemble(div_u1 * div_u1 * my_dx))))
    return u1

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
 def assemble_solution(self, t):  # returns
     """
     :param t: time
     :return: Womersley flow (analytic solution) at time t
     analytic solution at any time is a steady parabolic flow + linear combination of 8 modes
     modes were precomputed as 8 functions on given mesh and stored in hdf5 file
     """
     if self.tc is not None:
         self.tc.start('assembleSol')
     sol = Function(self.solutionSpace)
     # analytic solution has zero x and y components
     dofs2 = self.solutionSpace.sub(2).dofmap().dofs()  # gives field of indices corresponding to z axis
     sol.assign(Constant(("0.0", "0.0", "0.0")))  # QQ not needed
     sol.vector()[dofs2] += self.factor * self.bessel_parabolic.vector().array()  # parabolic part of sol
     for idx in range(8):  # add modes of Womersley sol
         sol.vector()[dofs2] += self.factor * cos(self.coefs_exp[idx] * pi * t) * self.bessel_real[idx].vector().array()
         sol.vector()[dofs2] += self.factor * -sin(self.coefs_exp[idx] * pi * t) * self.bessel_complex[idx].vector().array()
     if self.tc is not None:
         self.tc.end('assembleSol')
     return sol

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]

#.........这里部分代码省略.........
        info('  Initializing output files.')
        # for creating paraview scripts
        self.metadata['filename_base'] = self.problem_code + '_' + self.metadata['name']

        # assemble file dictionary
        if self.doSave:
            if self.doSaveDiff:
                self.fileDict.update(self.fileDictDiff)
            if self.metadata['hasTentativeV']:
                self.fileDict.update(self.fileDictTent)
                if self.doSaveDiff:
                    self.fileDict.update(self.fileDictTentDiff)
            if self.metadata['hasTentativeP']:
                self.fileDict.update(self.fileDictTentP)
                if self.doSaveDiff:
                    self.fileDict.update(self.fileDictTentPDiff)
        else:
            self.fileDict = {}
        if self.args.wss != 'none':
            self.fileDict.update(self.fileDictWSSnorm)
            if self.args.wss_method == 'expression':
                self.fileDict.update(self.fileDictWSS)
        if self.args.debug_rot:
            self.fileDict.update(self.fileDictDebugRot)
        # create and setup files
        for key, value in self.fileDict.iteritems():
            value['file'] = XDMFFile(mpi_comm_world(), self.str_dir_name + "/" + self.problem_code + '_' +
                                     self.metadata['name'] + value['name'] + ".xdmf")
            value['file'].parameters['rewrite_function_mesh'] = False  # save mesh only once per file

    # method for saving divergence (ensuring, that it will be one time line in ParaView)
    def save_div(self, is_tent, field):
        self.tc.start('div')
        self.divFunction.assign(project(div(field), self.divSpace))
        self.fileDict['d2' if is_tent else 'd']['file'].write(self.divFunction, self.actual_time)
        self.tc.end('div')

    def compute_div(self, is_tent, velocity):
        self.tc.start('divNorm')
        div_list = self.listDict['d2' if is_tent else 'd']['list']
        div_list.append(norm(velocity, 'Hdiv0'))
        self.tc.end('divNorm')

    # method for saving velocity (ensuring, that it will be one time line in ParaView)
    def save_vel(self, is_tent, field):
        self.vFunction.assign(field)
        self.fileDict['u2' if is_tent else 'u']['file'].write(self.vFunction, self.actual_time)
        if self.doSaveDiff:
            self.vFunction.assign((1.0 / self.vel_normalization_factor[0]) * (field - self.solution))
            self.fileDict['u2D' if is_tent else 'uD']['file'].write(self.vFunction, self.actual_time)

    def compute_err(self, is_tent, velocity, t):
        if self.doErrControl:
            er_list_L2 = self.listDict['u2L2' if is_tent else 'u_L2']['list']
            er_list_H1 = self.listDict['u2H1' if is_tent else 'u_H1']['list']
            self.tc.start('errorV')
            # assemble is faster than errornorm
            errorL2_sq = assemble(inner(velocity - self.solution, velocity - self.solution) * dx)
            errorH1seminorm_sq = assemble(inner(grad(velocity - self.solution), grad(velocity - self.solution)) * dx)
            info('  H1 seminorm error: %f' % sqrt(errorH1seminorm_sq))
            errorL2 = sqrt(errorL2_sq)
            errorH1 = sqrt(errorL2_sq + errorH1seminorm_sq)
            info("  Relative L2 error in velocity = %f" % (errorL2 / self.analytic_v_norm_L2))
            self.last_error = errorH1 / self.analytic_v_norm_H1
            self.last_status_functional = self.last_error
            info("  Relative H1 error in velocity = %f" % self.last_error)

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
    def solve(self, problem):
        self.problem = problem
        doSave = problem.doSave
        save_this_step = False
        onlyVel = problem.saveOnlyVel
        dt = self.metadata['dt']

        nu = Constant(self.problem.nu)
        # TODO check proper use of watches
        self.tc.init_watch('init', 'Initialization', True, count_to_percent=False)
        self.tc.init_watch('rhs', 'Assembled right hand side', True, count_to_percent=True)
        self.tc.init_watch('updateBC', 'Updated velocity BC', True, count_to_percent=True)
        self.tc.init_watch('applybc1', 'Applied velocity BC 1st step', True, count_to_percent=True)
        self.tc.init_watch('applybc3', 'Applied velocity BC 3rd step', True, count_to_percent=True)
        self.tc.init_watch('applybcP', 'Applied pressure BC or othogonalized rhs', True, count_to_percent=True)
        self.tc.init_watch('assembleMatrices', 'Initial matrix assembly', False, count_to_percent=True)
        self.tc.init_watch('solve 1', 'Running solver on 1st step', True, count_to_percent=True)
        self.tc.init_watch('solve 2', 'Running solver on 2nd step', True, count_to_percent=True)
        self.tc.init_watch('solve 3', 'Running solver on 3rd step', True, count_to_percent=True)
        self.tc.init_watch('solve 4', 'Running solver on 4th step', True, count_to_percent=True)
        self.tc.init_watch('assembleA1', 'Assembled A1 matrix (without stabiliz.)', True, count_to_percent=True)
        self.tc.init_watch('assembleA1stab', 'Assembled A1 stabilization', True, count_to_percent=True)
        self.tc.init_watch('next', 'Next step assignments', True, count_to_percent=True)
        self.tc.init_watch('saveVel', 'Saved velocity', True)

        self.tc.start('init')

        # Define function spaces (P2-P1)
        mesh = self.problem.mesh
        self.V = VectorFunctionSpace(mesh, "Lagrange", 2)  # velocity
        self.Q = FunctionSpace(mesh, "Lagrange", 1)  # pressure
        self.PS = FunctionSpace(mesh, "Lagrange", 2)  # partial solution (must be same order as V)
        self.D = FunctionSpace(mesh, "Lagrange", 1)   # velocity divergence space
        if self.bc == 'lagrange':
            L = FunctionSpace(mesh, "R", 0)
            QL = self.Q*L

        problem.initialize(self.V, self.Q, self.PS, self.D)

        # Define trial and test functions
        u = TrialFunction(self.V)
        v = TestFunction(self.V)
        if self.bc == 'lagrange':
            (pQL, rQL) = TrialFunction(QL)
            (qQL, lQL) = TestFunction(QL)
        else:
            p = TrialFunction(self.Q)
            q = TestFunction(self.Q)

        n = FacetNormal(mesh)
        I = Identity(u.geometric_dimension())

        # Initial conditions: u0 velocity at previous time step u1 velocity two time steps back p0 previous pressure
        [u1, u0, p0] = self.problem.get_initial_conditions([{'type': 'v', 'time': -dt},
                                                          {'type': 'v', 'time': 0.0},
                                                          {'type': 'p', 'time': 0.0}])

        if doSave:
            problem.save_vel(False, u0, 0.0)
            problem.save_vel(True, u0, 0.0)

        u_ = Function(self.V)         # current tentative velocity
        u_cor = Function(self.V)         # current corrected velocity
        if self.bc == 'lagrange':
            p_QL = Function(QL)    # current pressure or pressure help function from rotation scheme
            pQ = Function(self.Q)     # auxiliary function for conversion between QL.sub(0) and Q
        else:
            p_ = Function(self.Q)         # current pressure or pressure help function from rotation scheme
        p_mod = Function(self.Q)      # current modified pressure from rotation scheme

        # Define coefficients
        k = Constant(self.metadata['dt'])
        f = Constant((0, 0, 0))

        # Define forms
        # step 1: Tentative velocity, solve to u_
        u_ext = 1.5*u0 - 0.5*u1  # extrapolation for convection term

        # Stabilisation
        h = CellSize(mesh)
        # CBC delta:
        if self.cbcDelta:
            delta = Constant(self.stabCoef)*h/(sqrt(inner(u_ext, u_ext))+h)
        else:
            delta = Constant(self.stabCoef)*h**2/(2*nu*k + k*h*inner(u_ext, u_ext)+h**2)

        if self.use_full_SUPG:
            v1 = v + delta*0.5*k*dot(grad(v), u_ext)
            parameters['form_compiler']['quadrature_degree'] = 6
        else:
            v1 = v

        def nonlinearity(function):
            if self.use_ema:
               return 2*inner(dot(sym(grad(function)), u_ext), v1) * dx + inner(div(function)*u_ext, v1) * dx
                # return 2*inner(dot(sym(grad(function)), u_ext), v) * dx + inner(div(u_ext)*function, v) * dx
                # QQ implement this way?
            else:
                return inner(dot(grad(function), u_ext), v1) * dx

#.........这里部分代码省略.........

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]

#.........这里部分代码省略.........
        cell=cell_type
        )

    # Compute the problem
    errors = {'u': numpy.empty((len(mesh_sizes), len(Dt))),
              'p': numpy.empty((len(mesh_sizes), len(Dt)))
              }
    for k, mesh_size in enumerate(mesh_sizes):
        info('')
        info('')
        with Message('Computing for mesh size %r...' % mesh_size):
            mesh = mesh_generator(mesh_size)
            mesh_area = assemble(1.0 * dx(mesh))
            W = VectorFunctionSpace(mesh, 'CG', 2)
            P = FunctionSpace(mesh, 'CG', 1)
            method = MethodClass(W, P,
                                 rho, mu,
                                 theta=1.0,
                                 #theta=0.5,
                                 stabilization=None
                                 #stabilization='SUPG'
                                 )
            u1 = Function(W)
            p1 = Function(P)
            err_p = Function(P)
            divu1 = Function(P)
            for j, dt in enumerate(Dt):
                # Prepare previous states for multistepping.
                u = [Expression(
                    sol_u.cppcode,
                    degree=_truncate_degree(solution['u']['degree']),
                    t=0.0,
                    cell=cell_type
                    ),
                    # Expression(
                    #sol_u.cppcode,
                    #degree=_truncate_degree(solution['u']['degree']),
                    #t=0.5*dt,
                    #cell=cell_type
                    #)
                    ]
                sol_u.t = dt
                u_bcs = [DirichletBC(W, sol_u, 'on_boundary')]
                sol_p.t = dt
                #p_bcs = [DirichletBC(P, sol_p, 'on_boundary')]
                p_bcs = []
                fenics_rhs0.t = 0.0
                fenics_rhs1.t = dt
                method.step(dt,
                            u1, p1,
                            u, p0,
                            u_bcs=u_bcs, p_bcs=p_bcs,
                            f0=fenics_rhs0, f1=fenics_rhs1,
                            verbose=False,
                            tol=1.0e-10
                            )
                sol_u.t = dt
                sol_p.t = dt
                errors['u'][k][j] = errornorm(sol_u, u1)
                # The pressure is only determined up to a constant which makes
                # it a bit harder to define what the error is. For our
                # purposes, choose an alpha_0\in\R such that
                #
                #    alpha0 = argmin ||e - alpha||^2
                #
                # with  e := sol_p - p.
                # This alpha0 is unique and explicitly given by
                #
                #     alpha0 = 1/(2|Omega|) \int (e + e*)
                #            = 1/|Omega| \int Re(e),
                #
                # i.e., the mean error in \Omega.
                alpha = assemble(sol_p * dx(mesh)) \
                    - assemble(p1 * dx(mesh))
                alpha /= mesh_area
                # We would like to perform
                #     p1 += alpha.
                # To avoid creating a temporary function every time, assume
                # that p1 lives in a function space where the coefficients
                # represent actual function values. This is true for CG
                # elements, for example. In that case, we can just add any
                # number to the vector of p1.
                p1.vector()[:] += alpha
                errors['p'][k][j] = errornorm(sol_p, p1)

                show_plots = False
                if show_plots:
                    plot(p1, title='p1', mesh=mesh)
                    plot(sol_p, title='sol_p', mesh=mesh)
                    err_p.vector()[:] = p1.vector()
                    sol_interp = interpolate(sol_p, P)
                    err_p.vector()[:] -= sol_interp.vector()
                    #plot(sol_p - p1, title='p1 - sol_p', mesh=mesh)
                    plot(err_p, title='p1 - sol_p', mesh=mesh)
                    #r = Expression('x[0]', degree=1, cell=triangle)
                    #divu1 = 1 / r * (r * u1[0]).dx(0) + u1[1].dx(1)
                    divu1.assign(project(u1[0].dx(0) + u1[1].dx(1), P))
                    plot(divu1, title='div(u1)')
                    interactive()
    return errors

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
    def step(self,
             dt,
             u1, p1,
             u, p0,
             u_bcs, p_bcs,
             f0=None, f1=None,
             verbose=True,
             tol=1.0e-10
             ):
        # Some initial sanity checkups.
        assert dt > 0.0
        # Define trial and test functions
        ui = Function(self.W)
        v = TestFunction(self.W)
        # Create functions
        # Define coefficients
        k = Constant(dt)

        # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        # Compute tentative velocity step:
        #
        #     F(u) = 0,
        #     F(u) := rho (U0 + (u.\nabla)u) - mu \div(\nabla u) - f = 0.
        #
        with Message('Computing tentative velocity'):
            # TODO higher-order scheme for time integration
            #
            # For higher-order schemes, see
            #
            #     A comparison of time-discretization/linearization approaches
            #     for the incompressible Navier-Stokes equations;
            #     Volker John, Gunar Matthies, Joachim Rang;
            #     Comput. Methods Appl. Mech. Engrg. 195 (2006) 5995-6010;
            #     <http://www.wias-berlin.de/people/john/ELECTRONIC_PAPERS/JMR06.CMAME.pdf>.
            #
            F1 = self.rho * inner((ui - u[0])/k, v) * dx

            if abs(self.theta) > DOLFIN_EPS:
                # Implicit terms.
                if f1 is None:
                    raise RuntimeError('Implicit schemes need right-hand side '
                                       'at target step (f1).')
                F1 -= self.theta * _rhs_weak(ui, v, f1, self.rho, self.mu)
            if abs(1.0 - self.theta) > DOLFIN_EPS:
                # Explicit terms.
                if f0 is None:
                    raise RuntimeError('Explicit schemes need right-hand side '
                                       'at current step (f0).')
                F1 -= (1.0 - self.theta) \
                    * _rhs_weak(u[0], v, f0, self.rho, self.mu)

            if p0:
                F1 += dot(grad(p0), v) * dx

            #if stabilization:
            #    tau = stab.supg2(V.mesh(),
            #                     u_1,
            #                     mu/rho,
            #                     V.ufl_element().degree()
            #                     )
            #    R = rho*(ui - u_1)/k
            #    if abs(theta) > DOLFIN_EPS:
            #        R -= theta * _rhs_strong(ui, f1, rho, mu)
            #    if abs(1.0-theta) > DOLFIN_EPS:
            #        R -= (1.0-theta) * _rhs_strong(u_1, f0, rho, mu)
            #    if p_1:
            #        R += grad(p_1)
            #    # TODO use u_1 or ui here?
            #    F1 += tau * dot(R, grad(v)*u_1) * dx

            # Get linearization and solve nonlinear system.
            # If the scheme is fully explicit (theta=0.0), then the system is
            # actually linear and only one Newton iteration is performed.
            J = derivative(F1, ui)

            # What is a good initial guess for the Newton solve?
            # Three choices come to mind:
            #
            #    (1) the previous solution u_1,
            #    (2) the intermediate solution from the previous step ui_1,
            #    (3) the solution of the semilinear system
            #        (u.\nabla(u) -> u_1.\nabla(u)).
            #
            # Numerical experiments with the Karman vortex street show that the
            # order of accuracy is (1), (3), (2). Typical norms would look like
            #
            #     ||u - u_1 || = 1.726432e-02
            #     ||u - ui_1|| = 2.720805e+00
            #     ||u - u_e || = 5.921522e-02
            #
            # Hence, use u_1 as initial guess.
            ui.assign(u[0])

            # Wrap the solution in a try-catch block to make sure we call end()
            # if necessary.
            #problem = NonlinearVariationalProblem(F1, ui, u_bcs, J)
            #solver  = NonlinearVariationalSolver(problem)
            solve(
                F1 == 0, ui,
                bcs=u_bcs,
#.........这里部分代码省略.........

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
def solve_fixed_point(
        mesh,
        W_element, P_element, Q_element,
        u0, p0, theta0,
        kappa, rho, mu, cp,
        g, extra_force,
        heat_source,
        u_bcs, p_bcs,
        theta_dirichlet_bcs,
        theta_neumann_bcs,
        my_dx, my_ds,
        max_iter,
        tol
        ):
    # Solve the coupled heat-Stokes equation approximately. Do this
    # iteratively by solving the heat equation, then solving Stokes with the
    # updated heat, the heat equation with the updated velocity and so forth
    # until the change is 'small'.
    WP = FunctionSpace(mesh, MixedElement([W_element, P_element]))
    Q = FunctionSpace(mesh, Q_element)
    # Initialize functions.
    up0 = Function(WP)
    u0, p0 = up0.split()

    theta1 = Function(Q)
    for _ in range(max_iter):
        heat_problem = heat.Heat(
            Q,
            kappa=kappa,
            rho=rho(theta0),
            cp=cp,
            convection=u0,
            source=heat_source,
            dirichlet_bcs=theta_dirichlet_bcs,
            neumann_bcs=theta_neumann_bcs,
            my_dx=my_dx,
            my_ds=my_ds
            )

        theta1.assign(heat_problem.solve_stationary())

        # Solve problem for velocity, pressure.
        f = rho(theta0) * g  # coupling
        if extra_force:
            f += as_vector((extra_force[0], extra_force[1], 0.0))
        # up1 = up0.copy()
        stokes.stokes_solve(
            up0,
            mu,
            u_bcs, p_bcs,
            f,
            my_dx=my_dx,
            tol=1.0e-10,
            verbose=False,
            maxiter=1000
            )

        # from dolfin import plot
        # plot(u0)
        # plot(theta0)

        theta_diff = errornorm(theta0, theta1)
        info('||theta - theta0|| = {:e}'.format(theta_diff))
        # info('||u - u0||         = {:e}'.format(u_diff))
        # info('||p - p0||         = {:e}'.format(p_diff))
        # diff = theta_diff + u_diff + p_diff
        diff = theta_diff
        info('sum = {:e}'.format(diff))

        # # Show the iterates.
        # plot(theta0, title='theta0')
        # plot(u0, title='u0')
        # interactive()
        # #exit()
        if diff < tol:
            break

        theta0.assign(theta1)

    # Create a *deep* copy of u0, p0, to be able to deal with them as
    # actually separate entities.
    u0, p0 = up0.split(deepcopy=True)
    return u0, p0, theta0

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
    def step(self,
             dt,
             u1, p1,
             u, p0,
             u_bcs, p_bcs,
             f0=None, f1=None,
             verbose=True,
             tol=1.0e-10
             ):
        '''General pressure projection scheme as described in section 3.4 of
        :cite:`GMS06`.
        '''
        # Some initial sanity checkups.
        assert dt > 0.0

        # Define coefficients
        k = Constant(dt)

        # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        # Compute tentative velocity step.
        # rho (u[-1] + (u.\nabla)u) = mu 1/r \div(r \nabla u) + rho g.
        with Message('Computing tentative velocity'):
            solver = NewtonSolver()
            solver.parameters['maximum_iterations'] = 5
            solver.parameters['absolute_tolerance'] = tol
            solver.parameters['relative_tolerance'] = 0.0
            solver.parameters['report'] = True
            # The nonlinear term makes the problem generally nonsymmetric.
            solver.parameters['linear_solver'] = 'gmres'
            # If the nonsymmetry is too strong, e.g., if u[-1] is large, then
            # AMG preconditioning might not work very well.
            # Use HYPRE-Euclid instead of ILU for parallel computation.
            solver.parameters['preconditioner'] = 'hypre_euclid'
            solver.parameters['krylov_solver']['relative_tolerance'] = tol
            solver.parameters['krylov_solver']['absolute_tolerance'] = 0.0
            solver.parameters['krylov_solver']['maximum_iterations'] = 1000
            solver.parameters['krylov_solver']['monitor_convergence'] = verbose

            ui = Function(self.W)
            step_problem = \
                TentativeVelocityProblem(ui,
                                         self.theta,
                                         self.rho, self.mu,
                                         u, p0, k,
                                         u_bcs,
                                         f0, f1,
                                         stabilization=self.stabilization,
                                         dx=self.dx
                                         )

            # Take u[-1] as initial guess.
            ui.assign(u[-1])
            solver.solve(step_problem, ui.vector())
            #div_u = 1/r * div(r*ui)
            #plot(div_u)
            #interactive()
        # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        with Message('Computing pressure correction'):
            # The pressure correction is based on the update formula
            #
            #                           (    dphi/dr    )
            #     rho/dt (u_{n+1}-u*) + (    dphi/dz    ) = 0
            #                           (1/r dphi/dtheta)
            #
            # with
            #
            #     phi = p_{n+1} - p*
            #
            # and
            #
            #      1/r d/dr    (r ur_{n+1})
            #    +     d/dz    (  uz_{n+1})
            #    + 1/r d/dtheta(  utheta_{n+1}) = 0
            #
            # With the assumption that u does not change in the direction
            # theta, one derives
            #
            #  - 1/r * div(r * \nabla phi) = 1/r * rho/dt div(r*(u_{n+1} - u*))
            #  - 1/r * n. (r * \nabla phi) = 1/r * rho/dt  n.(r*(u_{n+1} - u*))
            #
            # In its weak form, this is
            #
            #   \int r * \grad(phi).\grad(q) *2*pi =
            #       - rho/dt \int div(r*u*) q *2*p
            #       - rho/dt \int_Gamma n.(r*(u_{n+1}-u*)) q *2*pi.
            #
            # (The terms 1/r cancel with the volume elements 2*pi*r.)
            # If Dirichlet boundary conditions are applied to both u* and u_n
            # (the latter in the final step), the boundary integral vanishes.
            #
            alpha = 1.0
            #alpha = 3.0 / 2.0
            rotational_form = False
            self._pressure_poisson(p1, p0,
                                   self.mu, ui,
                                   self.rho * alpha * ui / k,
                                   p_bcs=p_bcs,
                                   rotational_form=rotational_form,
                                   tol=tol,
                                   verbose=verbose
#.........这里部分代码省略.........

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]

#.........这里部分代码省略.........
                    f = interpolate(womersleyBC.analytic_pressure(self.factor, d['time']), self.pSpace)
            out.append(f)
        return out

    def get_outflow_measures(self):
        return [self.dsOut]

    def get_outflow_measure_form(self):
        return self.dsOut

    def get_v_solution(self, t):
        v = self.assemble_solution(t)
        return v

    def get_p_solution(self, t):
        p = interpolate(womersleyBC.analytic_pressure(self.factor, t), self.pSpace)
        return p

    def update_time(self, actual_time, step_number):
        super(Problem, self).update_time(actual_time, step_number)
        if self.actual_time > 0.5 and int(round(self.actual_time * 1000)) % 1000 == 0:
            self.isWholeSecond = True
            seconds = int(round(self.actual_time))
            self.second_list.append(seconds)
            self.N1 = seconds*self.stepsInSecond
            self.N0 = (seconds-1)*self.stepsInSecond
        else:
            self.isWholeSecond = False

        self.solution = self.assemble_solution(self.actual_time)

        # Update boundary condition
        self.tc.start('updateBC')
        self.v_in.assign(self.solution)
        self.tc.end('updateBC')

        # construct analytic pressure (used for computing pressure and force errors)
        self.tc.start('analyticP')
        analytic_pressure = womersleyBC.analytic_pressure(self.factor, self.actual_time)
        self.sol_p = interpolate(analytic_pressure, self.pSpace)
        self.tc.end('analyticP')

        self.tc.start('analyticVnorms')
        self.analytic_v_norm_L2 = norm(self.solution, norm_type='L2')
        self.analytic_v_norm_H1 = norm(self.solution, norm_type='H1')
        self.analytic_v_norm_H1w = sqrt(assemble((inner(grad(self.solution), grad(self.solution)) +
                                                  inner(self.solution, self.solution)) * self.dsWall))
        self.listDict['av_norm_L2']['list'].append(self.analytic_v_norm_L2)
        self.listDict['av_norm_H1']['list'].append(self.analytic_v_norm_H1)
        self.listDict['av_norm_H1w']['list'].append(self.analytic_v_norm_H1w)
        self.tc.end('analyticVnorms')

    def assemble_solution(self, t):  # returns Womersley sol for time t
        if self.tc is not None:
            self.tc.start('assembleSol')
        sol = Function(self.solutionSpace)
        dofs2 = self.solutionSpace.sub(2).dofmap().dofs()  # gives field of indices corresponding to z axis
        sol.assign(Constant(("0.0", "0.0", "0.0")))  # QQ not needed
        sol.vector()[dofs2] += self.factor * self.bessel_parabolic.vector().array()  # parabolic part of sol
        for idx in range(8):  # add modes of Womersley sol
            sol.vector()[dofs2] += self.factor * cos(self.coefs_exp[idx] * pi * t) * self.bessel_real[idx].vector().array()
            sol.vector()[dofs2] += self.factor * -sin(self.coefs_exp[idx] * pi * t) * self.bessel_complex[idx].vector().array()
        if self.tc is not None:
            self.tc.end('assembleSol')
        return sol

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
def test_time_step():
    problem = problems.Crucible()

    boundaries = problem.wp_boundaries

    # The melting point of GaAs is 1511 K.
    average_temp = 1520.0

    f = Constant(0.0)

    material = problem.subdomain_materials[problem.wpi]
    rho = material.density(average_temp)
    cp = material.specific_heat_capacity
    kappa = material.thermal_conductivity

    my_ds = Measure("ds")(subdomain_data=boundaries)

    # from dolfin import DirichletBC
    convection = None
    heat = maelstrom.heat.Heat(
        problem.Q,
        kappa,
        rho,
        cp,
        convection,
        source=Constant(0.0),
        dirichlet_bcs=problem.theta_bcs_d,
        neumann_bcs=problem.theta_bcs_n,
        robin_bcs=problem.theta_bcs_r,
        my_dx=dx,
        my_ds=my_ds,
    )

    # create time stepper
    # stepper = parabolic.ExplicitEuler(heat)
    stepper = parabolic.ImplicitEuler(heat)
    # stepper = parabolic.Trapezoidal(heat)

    theta0 = project(Constant(average_temp), problem.Q)
    # theta0 = heat.solve_stationary()
    theta0.rename("theta0", "temperature")

    theta1 = Function(problem.Q)
    theta1 = Function(problem.Q)

    t = 0.0
    dt = 1.0e-3
    end_time = 10 * dt

    with XDMFFile("temperature.xdmf") as f:
        f.parameters["flush_output"] = True
        f.parameters["rewrite_function_mesh"] = False

        f.write(theta0, t)
        while t < end_time:
            theta1.assign(stepper.step(theta0, t, dt))
            theta0.assign(theta1)
            t += dt
            #
            f.write(theta0, t)

    assert abs(maelstrom.helpers.average(theta0) - 1519.81) < 1.0e-2

    return

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
def compute_tentative_velocity(
    time_step_method, rho, mu, u, p0, dt, u_bcs, f, W, my_dx, tol
):
    """Compute the tentative velocity via

    .. math::
        \\rho (u_0 + (u\\cdot\\nabla)u) =
            \\mu \\frac{1}{r} \\div(r \\nabla u) + \\rho g.
    """

    class TentativeVelocityProblem(NonlinearProblem):
        def __init__(self, ui, time_step_method, rho, mu, u, p0, dt, bcs, f, my_dx):
            super(TentativeVelocityProblem, self).__init__()

            W = ui.function_space()
            v = TestFunction(W)

            self.bcs = bcs

            r = SpatialCoordinate(ui.function_space().mesh())[0]

            def me(uu, ff):
                return _momentum_equation(uu, v, p0, ff, rho, mu, my_dx)

            self.F0 = rho * dot(ui - u[0], v) / dt * 2 * pi * r * my_dx
            if time_step_method == "forward euler":
                self.F0 += me(u[0], f[0])
            elif time_step_method == "backward euler":
                self.F0 += me(ui, f[1])
            else:
                assert (
                    time_step_method == "crank-nicolson"
                ), "Unknown time stepper '{}'".format(
                    time_step_method
                )
                self.F0 += 0.5 * (me(u[0], f[0]) + me(ui, f[1]))

            self.jacobian = derivative(self.F0, ui)
            self.reset_sparsity = True
            return

        # pylint: disable=unused-argument
        def F(self, b, x):
            # We need to evaluate F at x, so we have to make sure that self.F0
            # is assembled for ui=x. We could use a self.ui and set
            #
            #     self.ui.vector()[:] = x
            #
            # here. One way around this copy is to instantiate this class with
            # the same Function ui that is then used for the solver.solve().
            assemble(self.F0, tensor=b, form_compiler_parameters={"optimize": True})
            for bc in self.bcs:
                bc.apply(b, x)
            return

        def J(self, A, x):
            # We can ignore x; see comment at F().
            assemble(
                self.jacobian, tensor=A, form_compiler_parameters={"optimize": True}
            )
            for bc in self.bcs:
                bc.apply(A)
            self.reset_sparsity = False
            return

    solver = NewtonSolver()
    solver.parameters["maximum_iterations"] = 10
    solver.parameters["absolute_tolerance"] = tol
    solver.parameters["relative_tolerance"] = 0.0
    solver.parameters["report"] = True
    # While GMRES+ILU converges if the time step is small enough, increasing
    # the time step slows down convergence dramatically in some cases. This
    # makes the step fail, and the adaptive time stepper will decrease the step
    # size. This size can be _very_ small such that simulation take forever.
    # For now, just use a direct solver. Choose UMFPACK over SuperLU since the
    # docker image doesn't contain SuperLU yet, cf.
    # <https://bitbucket.org/fenics-project/docker/issues/64/add-superlu>.
    # TODO come up with an appropriate GMRES preconditioner here
    solver.parameters["linear_solver"] = "umfpack"

    ui = Function(W)
    step_problem = TentativeVelocityProblem(
        ui, time_step_method, rho, mu, u, p0, dt, u_bcs, f, my_dx
    )

    # Take u[0] as initial guess.
    ui.assign(u[0])
    solver.solve(step_problem, ui.vector())

    # Make sure ui is from W. This should happen anyways, but somehow doesn't.
    # TODO find out why not
    ui = project(ui, W)
    # div_u = 1/r * div(r*ui)
    return ui

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import assign [as 别名]
class OSI(Field):

    @classmethod
    def default_params(cls):
        params = Field.default_params()
        params.update(
            finalize=True,
            )
        return params

    def add_fields(self):
        params = self.params.copy_recursive()
        params["save"] = False
        params["plot"] = False
        #params["callback"] = False
        #params.pop("finalize")

        fields = []
        #fields.append(WSS(params=params))
        #return fields

        f = TimeIntegral("WSS", params=params, label="OSI")
        fields.append(f)
        fields.append(Magnitude(f, params=params))

        f = Magnitude("WSS", params=params)
        fields.append(f)
        fields.append(TimeIntegral(f, params=params, label="OSI"))

        #f = TimeIntegral("WSS", label="OSI")
        #fields.append(f)
        #fields.append(Magnitude(f))

        #f = Magnitude("WSS")
        #fields.append(f)
        #fields.append(TimeIntegral(f, label="OSI"))

        return fields

    def before_first_compute(self, get):
        tau = get("WSS")
        self.osi = Function(tau.sub(0).function_space().collapse())

    def compute(self, get):
        # Requires the fields Magnitude(TimeIntegral("WSS", label="OSI")) and
        # TimeIntegral(Magnitude("WSS"), label="OSI")
        #self.mag_ta_wss = get("Magnitude_TimeIntegral_WSS_OSI")
        #self.ta_mag_wss = get("TimeIntegral_Magnitude_WSS_OSI")
        self.mag_ta_wss = get("Magnitude_TimeIntegral_WSS-OSI")
        self.ta_mag_wss = get("TimeIntegral_Magnitude_WSS-OSI")

        if self.params.finalize:
            return None
        elif self.mag_ta_wss == None or self.ta_mag_wss == None:
            return None
        else:
            expr = conditional(self.ta_mag_wss < 1e-15,
                               0.0,
                               0.5 * (1.0 - self.mag_ta_wss / self.ta_mag_wss))
            self.osi.assign(project(expr, self.osi.function_space()))
            return self.osi

    def after_last_compute(self, get):
        self.mag_ta_wss = get("Magnitude_TimeIntegral_WSS-OSI")
        self.ta_mag_wss = get("TimeIntegral_Magnitude_WSS-OSI")
        #print self.name, " Calling after_last_compute"

        expr = conditional(self.ta_mag_wss < 1e-15,
                           0.0,
                           0.5 * (1.0 - self.mag_ta_wss / self.ta_mag_wss))
        self.osi.assign(project(expr, self.osi.function_space()))

        return self.osi