Add a vector to a scalar - Forward Eular method

1

I apologize for my earlier untidy description. Here is a better version of my question

I’m trying to use Forward Eular method in solving the heat equation mentioned here

Then, the PDE becomes:

\frac{u^{n+1}-u^n}{\Delta t}=\triangledown^2u^n+f^n

And by following a similar calculation as given in the tutorial, we get:

\int\limits_{\Omega}uv dx=\int\limits_{\Omega}[u^n+\Delta t(\triangledown^2u^n+f^n)]v dx

Thus
a(u,v)=\int\limits_{\Omega}uv dx
L(v)=\int\limits_{\Omega}[u^n+\Delta t(\triangledown^2u^n+f^n)]v dx

where f is a constant function.

Now to compute the variational form in the code, I will have to compute the \triangledown^2u^n+f^n part which I’m finding it difficult.

Obviously the ufl.grad(ufl.grad(u_n))+f is not going to work because of the mismatch in the dimensions.
I first thought about having f as a constant numpy array. But I believe there should be a better way of handling this issue. And even for the numpy array I’m not exactly sure how to get the corresponding dimension related to ufl.grad(ufl.grad(u_n))

Appreciate your help and sorry again for the previous untidy question.

2

First of all, please format your code with markdown syntax, ie

```python
def f(x):
    return x
```

do you want to add f to each component of the matrix grad(grad(u_n))?

Please not that the following does not make sense, as double_derivativeis a 2x2 matrix, v a scalar value.

Please write down your variational form in latex formatting to make it easier for people to understand what you are trying to do.

3

I’m extremely sorry for the untidy format. I have edited it in a better way now

4

\Delta^2u is not equal to ufl.grad(ufl.grad(u)), But ufl.div(ufl.grad(u)).
Secondly, in the Finite element method, one uses integration by parts to reduce this term to `ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx. See for instance:
https://jorgensd.github.io/dolfinx-tutorial/chapter1/fundamentals.html#finite-element-variational-formulation

1 Like
5

Thank you so much. With your help I managed to do that.
So, L(v)=\int\limits_{\Omega}u^nvdx-\Delta t\int\limits_{\Omega}\triangledown u^n\triangledown v dx+\Delta t\int\limits_{\Omega}f^nvdx

I’m posting my entrie code below just in case if anyone else had the same problem:

import numpy
from dolfinx import mesh, fem
import ufl
from mpi4py import MPI
from petsc4py import PETSc

t = 0 # Start time
T = 2 # End time
num_steps = 20 # Number of time steps
dt = (T-t)/num_steps # Time step size
alpha = 3
beta = 1.2

nx, ny = 5, 5
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("CG", 1))

class exact_solution():
    def __init__(self, alpha, beta, t):
        self.alpha = alpha
        self.beta = beta
        self.t = t
    def __call__(self, x):
        return 1 + x[0]**2 + self.alpha * x[1]**2 + self.beta * self.t
u_exact = exact_solution(alpha, beta, t)

u_D = fem.Function(V)
u_D.interpolate(u_exact)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))

u_n = fem.Function(V)
u_n.interpolate(u_exact)

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

f = fem.Constant(domain, beta - 2 - 2 * alpha)

def solve_problem(a,L,V,u_exact,num_steps,u_D):
    A = fem.petsc.assemble_matrix(a, bcs=[bc])
    A.assemble()
    b = fem.petsc.create_vector(L)
    uh = fem.Function(V)

    solver = PETSc.KSP().create(domain.comm)
    solver.setOperators(A)
    solver.setType(PETSc.KSP.Type.PREONLY)
    solver.getPC().setType(PETSc.PC.Type.LU)


    for n in range(num_steps):
        # Update Diriclet boundary condition 
        u_exact.t+=dt
        u_D.interpolate(u_exact)

        # Update the right hand side reusing the initial vector
        with b.localForm() as loc_b:
            loc_b.set(0)
        fem.petsc.assemble_vector(b, L)

        # Apply Dirichlet boundary condition to the vector
        fem.petsc.apply_lifting(b, [a], [[bc]])
        b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        fem.petsc.set_bc(b, [bc])

        # Solve linear problem
        solver.solve(b, uh.vector)
        uh.x.scatter_forward()

        # Update solution at previous time step (u_n)
        u_n.x.array[:] = uh.x.array

        # Compute L2 error and error at nodes
    V_ex = fem.FunctionSpace(domain, ("CG", 2))
    u_ex = fem.Function(V_ex)
    u_ex.interpolate(u_exact)
    error_L2 = numpy.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((uh - u_ex)**2 * ufl.dx)), op=MPI.SUM))
    if domain.comm.rank == 0:
        print(f"L2-error: {error_L2:.2e}")

    # Compute values at mesh vertices
    error_max = domain.comm.allreduce(numpy.max(numpy.abs(uh.x.array-u_D.x.array)), op=MPI.MAX)
    if domain.comm.rank == 0:
        print(f"Error_max: {error_max:.2e}")


a_backward_Euler = fem.form(u * v * ufl.dx + dt*ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx)
L_backward_Euler = fem.form((u_n + dt * f) * v * ufl.dx)
solve_problem(a_backward_Euler,L_backward_Euler,V,u_exact,num_steps,u_D)


a_forward_Euler = fem.form(u * v * ufl.dx)
L_forward_Euler = fem.form((u_n * v * ufl.dx)-(dt*ufl.dot(ufl.grad(u_n),ufl.grad(v))*ufl.dx)+(dt*f*v*ufl.dx))
solve_problem(a_backward_Euler,L_backward_Euler,V,u_exact,num_steps,u_D)

Source: Heat equation -Analytic solution