KSP solver issue

dolfinx
1

I am having trouble running my heat diffusion simulation. PETSc error code 73 keeps appearing when running the simulation, saying that “[0] Object is in wrong state” and “[0] Not for unassembled matrix”, referring to the line which says:

solver.solve(b_vec, T_new)

Can someone please help? The full script code is given below.

import dolfinx
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.io import XDMFFile
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, locate_dofs_geometrical, Constant, form)
from dolfinx import mesh, fem
import pyvista as pv


def visualiseMesh(msh: mesh):
    '''
    Provides a 2D rendering file (.xdmf) to show the mesh structure.  
    '''
    print("Plotting mesh in 3D using PyVista...")

    # Create XDMF file to save the mesh
    with XDMFFile(msh.comm, "mesh.xdmf", "w") as xdmf:
        xdmf.write_mesh(msh)

    # Extract the mesh topology (cell connectivity) and geometry (vertex coordinates)
    num_cells = msh.topology.index_map(msh.topology.dim).size_local
    cell_connectivity = msh.topology.connectivity(2, 0).array

    # Create the VTK-style connectivity array where the first value for each cell is the number of vertices
    vtk_cells = np.hstack([[3, *cell_connectivity[i:i+3]] for i in range(0, len(cell_connectivity), 3)])

    # Set the VTK cell types (5 corresponds to triangles)
    cell_types = np.full(num_cells, 5, dtype=np.uint8)

    # Get the vertex coordinates (geometry)
    points = msh.geometry.x

    # Convert the mesh into PyVista's UnstructuredGrid format
    grid = pv.UnstructuredGrid(vtk_cells, cell_types, points)

    # Plot the mesh using PyVista
    plotter = pv.Plotter()
    plotter.add_mesh(grid, show_edges=True)
    plotter.show()


def findHeatSolution(length, width, rho_param, c_p_param, k_param ):
    print("""
        ===================
        Starting simulation
        ===================
          """)
    
    print("Generating mesh...\n")
    # Creating the mesh (assuming a rectangular substrate for now)
    msh = mesh.create_rectangle(
        comm=MPI.COMM_WORLD,
        points=((0.0,0.0), (length, width)),
        n = (32, 16),
        cell_type=mesh.CellType.triangle,
    )

    #visualiseMesh(msh)

    # Defining the function space V
    V = fem.FunctionSpace(msh, ("Lagrange", 1))

    # Defining trial and test functions u and v
    T_n = Function(V) # Temperature at previous time step
    v = ufl.TestFunction(V)
    u = ufl.TrialFunction(V) # u is the temperature at the current timestep (what we are trying to solve)
    n = ufl.FacetNormal(msh)

    delta_t = 0.1

    # Defining the localised heat source power dissipation (uW) and coordinates
    # STILL NEED TO CHANGE THESE VALUES TO BE READ IN, RATHER THAN HARDCODED
    # print("Defining the localised heat source power dissipation and coordinates... \n")
    P1, P2, P3 = 420, 250, 500

    x1, y1 = -50.0, 0.0
    x2, y2 = 0.0, 0.0
    x3, y3 = 50.0, 0.0

    def localized_heat_source(x, power1=P1, power2=P2, power3=P3, radius=0.05):
        return (
        power1 * np.exp(-((x[0] - x1) ** 2 + (x[1] - y1) ** 2) / (2 * radius ** 2))
        + power2 * np.exp(-((x[0] - x2) ** 2 + (x[1] - y2) ** 2) / (2 * radius ** 2))
        + power3 * np.exp(-((x[0] - x3) ** 2 + (x[1] - y3) ** 2) / (2 * radius ** 2))
        )
    
   # Define the heat source as a Function in the FunctionSpace
    Q = Function(V)
    Q.interpolate(localized_heat_source)

    # Define the bilinear form (should contain unknowns u and v)
    a = (rho_param * c_p_param / delta_t) * u * v * ufl.dx + k_param * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - k_param * ufl.dot(ufl.grad(u), n) * v * ufl.ds

    # Define the linear form (including contributions from known T_n and Q)
    L = (rho_param * c_p_param / delta_t) * T_n * v * ufl.dx + Q * v * ufl.dx

    # Apply boundary conditions (e.g., Dirichlet with T=50 on all boundaries)
    boundary_conditions = [dirichletbc(PETSc.ScalarType(50.0), locate_dofs_geometrical(V, lambda x: np.full(x.shape[1], True)), V)]

    # Time-stepping parameters
    num_steps = 50  # Number of time steps
    time = 0.0

    with XDMFFile(MPI.COMM_WORLD, "heat_solution.xdmf", "w") as file:
        file.write_mesh(msh)
        file.write_function(T_n, time)

        # Assemble the matrix A
        A = fem.petsc.create_matrix(form(a))
        A_mat = fem.petsc.assemble_matrix(A, form(a), bcs=boundary_conditions)

        # Create a solver
        solver = PETSc.KSP().create(msh.comm)
        solver.setOperators(A_mat)  # Use the PETSc matrix
        solver.setType(PETSc.KSP.Type.CG)
        solver.getPC().setType(PETSc.PC.Type.LU)
        solver.setFromOptions()  # This is necessary to finalize the setup of the solver

        for step in range(num_steps):
            time += delta_t
            print(f"Time step {step + 1}/{num_steps}, Time: {time:.2f}")

            # Assemble the right-hand side vector (b)
            b = fem.petsc.create_vector(form(L))
            b_vec = fem.petsc.assemble_vector(b, form(L))
            
            # Apply boundary conditions to the vector b_vec
            fem.apply_lifting(b_vec, [form(a)], bcs=[boundary_conditions])
            b_vec.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
            fem.set_bc(b_vec, boundary_conditions)

            # Create the solution vector
            T_new = Function(V).vector
            
            #print("\nSize of b_vec:", b_vec.size)
            #print("Size of T_new.vector:", T_new.vector.size)

            solver.solve(b_vec, T_new)
            T_new.x.scatter_forward()

            # Update T_n for the next time step
            T_n.x.array[:] = T_new.x.array[:]

            # Save the current solution
            file.write_function(T_new, time)


    print("""
        ===================
        Simulation complete
        ===================
          """)
2

Add A_mat.assemble() (or A.assemble() as it is the same) after this call.

1 Like