How to do a convolution

117 views
Skip to first unread message

nathan magnan

unread,
Oct 10, 2024, 4:50:08 AM10/10/24
to Dedalus Users
Hello,

I am trying to solve an EVP, but one of the internal equations requires a convolution. Mathematically, it looks like this:
variable_2 (x) = integral[ dx' Kernel(x, x') variable_1(x') ]

This convolution can be conceived as an operator acting on variable_1, but I don't know how to write it in Dedalus.

As a side note, if omega is my eigenvalue, is omega^2 allowed to appear in one of the internal equations, or would that break linearity?

Thanks for your help !
Nathan

nathan magnan

unread,
Oct 10, 2024, 10:47:40 AM10/10/24
to Dedalus Users
Hello again,

I've tried to build a work-around that looks like this:

# i) build the trial functions as fields
Trial_functions = []

for i in range(Nx):
    name = 'trial_function_' + str(i)
    trial_function = dist.Field(name = name, bases = basis)
    trial_function['c'][i] = 1
   
    Trial_functions.append(trial_function)


# ii) build x = cst slices of the kernel as fields

def kernel(x_outer, x_inner):
    return # WHATEVER MY KERNEL IS

K_slices = []

for i in range(Nx):
    x_outer = x[i]
    name = 'K_slice_' + str(i)
    K_slice = dist.Field(name = name, bases = basis)

    for j in range(Nr):
        x_inner = r[j]
        K_slice['g'][j] = kernel(x_outer, x_inner)

    K_slices.append(K_slice)


# iii) build the image of the trial functions via my convolution operator. The output should be fields

Images_of_trial_functions = []

for i in range(Nx):
    trial_function = Trial_functions[i]
   
    name = 'image_of_trial_function_' + str(i)
    image_of_trial_function = dist.Field(name = name, bases = basis)

    for j in range(Nx):
        K_slice = K_slices[j]

        image_of_trial_function['g'][j] = (d3.Integrate(K_slice * trial_function, 'x').evaluate())['g'][0]
        K_slice.change_scales(1)
        trial_function.change_scales(1)
        image_of_trial_function.change_scales(1)
   
    Images_of_trial_functions.append(image_of_trial_function)


## iv) build the operator

variable_2 = 0 * variable_1

for i in range(Nx):
    trial_function = Trial_functions[i]
    image_of_trial_function = Images_of_trial_functions[i]

    coefficient_of_Sigma = d3.Integrate(trial_function * variable_1, 'x') # I DON'T KNOW IF THIS IS EQUIVALENT TO variable_1['c'][i]
    variable_2 += coefficient_of_Sigma * image_of_trial_function

My first issue concerns the penultimate line.
I'm not sure that the i-th element in variable_1['c'] is always given by integral(variable_1 * trial_function_i).
I'm afraid that this only works for some choices of trial functions. For instance, the trial function may form a basis that is not orthonormal, or they might be orthonormal but the scalar product may not be the direct integral, there could be a complex conjugation hidden in there.
Ideally, I would like an operator that gives me the i-th coefficient in the decomposition of variable_1. A bit like d3.Interpolate gives me the k-th coefficient in the grid representation of variable_1. Is such a thing available ? If not, any idea for a workaround?

My second issue is that I'm worried about the numerical efficiency of my method. How catastrophicly slow is it gonna be?

Thanks for your help !
Nathan

Daniel Lecoanet

unread,
Oct 10, 2024, 11:46:27 AM10/10/24
to Dedalus Users
Hi Nathan,

If you’re doing a periodic problem, you may be able to use that a convolution in physical space is a product in wavenumber space. Then you would want to make a GeneralFunction operator which multiplies by the Fourier transform of your kernel in coefficient space.

I haven’t thought about convolutions in bounded domains.

Daniel

-- 
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/f2cb215d-8ea2-4746-93e8-97f4f82d39e0n%40googlegroups.com.

Nathan Magnan

unread,
Oct 10, 2024, 11:48:22 AM10/10/24
to dedalu...@googlegroups.com
Hi Daniel,

Thanks for the idea!
Unfortunately, I have Dirichlet boundary conditions on both sides.

Best regards,
Nathan
You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/kEhMmyv2-Dw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/340729F4-C381-4AB2-87D7-67F213C13593%40northwestern.edu.

Reply all
Reply to author
Forward
0 new messages