Setting values to Neumann boundary

General
1

Hi everyone. I’m trying to solve the Poisson problem in a spherical shell; the outer sphere is a Dirichlet boundary, and the inner sphere is a Neumann boundary.
I want to specify the du/dn distribution at the inner boundary. I tried spherical harmonics to construct the distribution. The function is (I use sympy because I cannot find legendre in ufl package):

def magmap_IB(x,sph_file='SPHs/mrmqc_c2074.dat'):
    # read the spherical harmonic coefficients
    sphs_df = pd.read_csv(sph_file, header=None, sep='\s+', names=['l', 'm', 'glm', 'hlm'])
    glm_array = np.zeros((61, 61))
    hlm_array = np.zeros((61, 61))
    for i in range(sphs_df.shape[0]):
        glm_array[sphs_df['l'][i], sphs_df['m'][i]] = sphs_df['glm'][i]
        hlm_array[sphs_df['l'][i], sphs_df['m'][i]] = sphs_df['hlm'][I]
    
    import sympy as sp
    cos_colat = x[2] / (x[0] ** 2 + x[1] ** 2 + x[2] ** 2) ** (1 / 2)
    lon = ufl.acos(x[0] / ufl.sqrt(x[0] ** 2 + x[1] ** 2)) * ufl.sign(x[1]) + ufl.pi - ufl.sign(x[1]) * ufl.pi

    l_max = 17

    COSCOLAT,LON = sp.symbols('cos_colat,lon')
    Br = 0
    from ufl import sqrt,sin,cos
    Br = sp.symbols('0.')
    for i_l in range(l_max + 1):
        for i_m in range(i_l + 1):
            print('l: ', i_l, '; m:', i_m)
            Plm = sp.assoc_legendre(i_l, i_m, COSCOLAT)\
                  * (2 * math.factorial(i_l - i_m)
                     / math.factorial(i_l + i_m)) ** (1 / 2)

            Br = Br + Plm * (glm_array[i_l, i_m] * sp.cos(i_m * LON)
                             + hlm_array[i_l, i_m] * sp.sin(i_m * LON))
    magmap = eval(str(Br))
    return magmap

The main code is:

    V = fem.FunctionSpace(msh, ('CG', 3))
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    x = ufl.SpatialCoordinate(msh)

    # 边界自由度
    inner_boundary_dof = fem.locate_dofs_topological(V=V, entity_dim=2, entities=inner_boundary)
    outer_boundary_dof = fem.locate_dofs_topological(V=V, entity_dim=2, entities=outer_boundary)

    phi0 = 0.
    B0 = 100. * 1e-9
    mu0 = 4 * np.pi * 1e-7

    f = fem.Constant(msh, ScalarType(0))
    # dirichlet B.C.
    bc = fem.dirichletbc(value=ScalarType(phi0), dofs=outer_boundary_dof, V=V)

    magmap = magmap_IB(x)
   
    a = dot(grad(u), grad(v)) * dx
    L = dot(f, v) * dx - dot(magmap, v) * ds

    # Solve
    ksp_type = 'gmres'
    problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={'ksp_type': ksp_type})
    uh = problem.solve()

My problems are:
(1) The code cannot work with spherical harmonics over L=17 (the error is as follows). Are there any functions in dolfinx that can generate and handle with legendre polynomials?

libffcx_forms_cb3037c64ed6c464f7c84f3d41225a03793c8870.c:59428:26: error: integer literal is too large to be represented in any integer type
    sv_e0a[223] = -1 * ((495179769008019818390136611716089140625 * sv_e0a[2]) / 2) + (9738535457157723095006020030416419765625 * pow(sv_e0a[2], 3)) / 2;

(2) My plan B is to directly interpolate the distribution to the inner boundary. Is there a way in dolfinx to do that?

2

As you haven’t provided the CSV that you are using or the mesh, it is quite hard to give you any guidance.

3

Thanks a lot for replying! I upload the files to google drive.

The mesh: https://drive.google.com/file/d/13u2Z8sCbG6Uzf1gmzJhVcaWCL4wSI28M/view?usp=drive_link

The data distribution for the inner boundary (360*180): https://drive.google.com/file/d/1aGFLtr5VzdkmBaeMlFDE7D-JrelwLzJP/view?usp=drive_link

The spherical harmonic coefficients (file used in function magmap_IB):
https://drive.google.com/file/d/1MHJx7mkLobY4QK8XzrxHrtg0LzHQFAs5/view?usp=drive_link

4

Could you please make your script complete, with all import statements and loading of the mesh.

5

Here’s my full code:

import numpy as np
import pyvista
import ufl
from dolfinx import fem, io, plot
from dolfinx.io.gmshio import read_from_msh
from ufl import ds, dx, grad, dot

import pandas as pd

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import math

def xyz2rlonlat(xyz, for_psi=False):
    r_carrington = (xyz[0] ** 2 + xyz[1] ** 2 + xyz[2] ** 2)**(1/2)

    lon_carrington = math.asin(xyz[1] / (xyz[0] ** 2 + xyz[1] ** 2)**(1/2))
    if xyz[0] < 0:
        lon_carrington = np.pi - lon_carrington
    if lon_carrington < 0:
        lon_carrington += 2 * np.pi

    lat_carrington = np.pi / 2 - math.acos(xyz[2] / r_carrington)
    if for_psi:
        lat_carrington = np.pi / 2 - lat_carrington
    return [r_carrington, lon_carrington, lat_carrington]


def magmap_IB(x,sph_file='SPHs/mrmqc_c2158.dat'):
    sphs_df = pd.read_csv(sph_file, header=None, sep='\s+', names=['l', 'm', 'glm', 'hlm'])
    glm_array = np.zeros((61, 61))
    hlm_array = np.zeros((61, 61))
    for i in range(sphs_df.shape[0]):
        glm_array[sphs_df['l'][i], sphs_df['m'][i]] = sphs_df['glm'][i]
        hlm_array[sphs_df['l'][i], sphs_df['m'][i]] = sphs_df['hlm'][i]
    import sympy as sp

    cos_colat = x[2] / (x[0] ** 2 + x[1] ** 2 + x[2] ** 2) ** (1 / 2)
    lon = ufl.acos(x[0] / ufl.sqrt(x[0] ** 2 + x[1] ** 2)) * ufl.sign(x[1]) + ufl.pi - ufl.sign(x[1]) * ufl.pi

    l_max = 60

    COSCOLAT,LON = sp.symbols('cos_colat,lon')
    Br = 0
    from ufl import sqrt,sin,cos
    Br = sp.symbols('0.')
    for i_l in range(l_max + 1):
        for i_m in range(i_l + 1):
            print('l: ', i_l, '; m:', i_m)
            Plm = sp.assoc_legendre(i_l, i_m, COSCOLAT)\
                  * (2 * math.factorial(i_l - i_m)
                     / math.factorial(i_l + i_m)) ** (1 / 2)

            Br = Br + Plm * (glm_array[i_l, i_m] * sp.cos(i_m * LON)
                             + hlm_array[i_l, i_m] * sp.sin(i_m * LON))
    magmap = eval(str(Br))
    return magmap


def fem_laplace_in_shell(path,filename,inner_boundary_marker,outer_boundary_marker):
    msh, cell_tags, facet_tags = read_from_msh(path + filename + '.msh', MPI.COMM_WORLD, 0, gdim=3)
    inner_boundary = facet_tags.find(inner_boundary_marker)
    outer_boundary = facet_tags.find(outer_boundary_marker)

    # %% 创建函数空间，
    V = fem.FunctionSpace(msh, ('CG', 3))
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    x = ufl.SpatialCoordinate(msh)

    inner_boundary_dof = fem.locate_dofs_topological(V=V, entity_dim=2, entities=inner_boundary)
    outer_boundary_dof = fem.locate_dofs_topological(V=V, entity_dim=2, entities=outer_boundary)

    phi0 = 0.
    B0 = 100. * 1e-9
    mu0 = 4 * np.pi * 1e-7
    Rss = 2.5

    f = fem.Constant(msh, ScalarType(0))

    bc = fem.dirichletbc(value=ScalarType(phi0), dofs=outer_boundary_dof, V=V)

    # dipole field
    # u_dipole = -B0 / mu0 * x[2] / (x[0] ** 2 + x[1] ** 2 + x[2] ** 2) ** (1 / 2)

    magmap = magmap_IB(x)

    a = dot(grad(u), grad(v)) * dx
    L = dot(f, v) * dx - dot(magmap, v) * ds

    # Solve
    ksp_type = 'gmres'
    from datetime import datetime
    d1 = datetime.now()
    problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={'ksp_type': ksp_type})
    # problem.solver.view()
    uh = problem.solve()
    print('Solved')
    d2 = datetime.now()
    print(d2 - d1)

    # SAVE
    # with io.XDMFFile(msh.comm, "out_potentialfield/second_dipole_" + ksp_type + '_' + filename + ".xdmf", "w") as file:
    #     file.write_mesh(msh)
    #     file.write_function(uh)

    cells, types, x = plot.create_vtk_mesh(V)
    result = pyvista.UnstructuredGrid(cells, types, x)
    result.point_data["u"] = uh.x.array.real
    return result


if __name__ == '__main__':
    # Read MESH
    path = 'MESH/PHYSICAL/'
    filename = '(SphO1-SphI1)phy'

    # Extract boundaries
    inner_boundary_marker = 2
    outer_boundary_marker = 1

    result = fem_laplace_in_shell(path,filename,inner_boundary_marker,outer_boundary_marker)
    # %% Visualization
    result.set_active_scalars("u")
    plotter = pyvista.Plotter()
    plotter.add_mesh(result, show_edges=True, opacity=0.5, clim=[-0.1, 0.1], cmap='seismic')
    plotter.add_axes()
    plotter.show_grid()
    plotter.show()

Thanks again!