Multiplying a sum of gaussian source terms by decision variables in a poisson-constrained pde optimisation problem

Question

Hello, I am new to using FEniCS and dolfin-adjoint and am trying to solve a continuous version of the source inversion problem in two dimensions. I am having quite a bit of difficulty trying to create the gaussian source terms, as I need to multiply each of K total source terms with a continuous deicison variable in the weak formulation of the poisson equation, after which I need to solve for the decision variables, (and only allowing S out of K of the variables to be nonzero). I have been unable to find any sources on how to formulate these source terms--my code is below, and I think everything else works, except for the weak formulation which I am having trouble creating. Any advice or help would be much appreciated!

from __future__ import print_function
from fenics import *
from dolfin import *
from dolfin_adjoint import *
import numpy as np
#import pyipopt as opt
import random
import moola


#Using an nxn square mesh and S total turned on source terms
S = 3
n = 20
mesh = UnitSquareMesh(n, n)

# Function Space

V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "CG", 1)

# Create gaussian source term function and target locations; choose a value for constant sigma
sigma = .2
def source_term(xk,yk, sigma, degree, Fspace):
    return interpolate(Expression('exp(-sqrt(pow(x[0] - xk,2) + pow(x[1] - yk,2)) / 
            pow(sigma, 2))', xk=xk, yk = yk, sigma=sigma, degree = degree), Fspace)

   # Dirichlet boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, "on_boundary")

# Choosing sources as grid (this is just for visualization purposes)
nSize = np.linspace(1,n -1,n-1)
sourceGrid = []
for a in nSize:
    for b in nSize:
        sourceGrid.append([a,b])



# A list of all source terms, one for each w_kl with source centered at x_kl, y_kl 
source_expression_list = [source_term(xk, yk, sigma, 2, W) for xk, yk in sourceGrid]

# Each term in list is f_kl * v * dx
u = Function(V)
v = TestFunction(V)

# I would like to multiply each term in the list by a decision variable
# I am not sure how to do this in terms of syntax
rhs_weakForm = [assemble(-f * v * dx) for f in source_expression_list]

#I need to multiply the RHS of the weak formulation by the decision variables
F = inner(grad(v), grad(u))*dx - rhs_weakForm
solve(F == 0, u, bc, solver_parameters={"newton_solver":
                   {"absolute_tolerance": 1.0e-7,    "maximum_iterations": 20}})


if __name__ == "__main__":
    x_init = np.zeros(len(sourceGrid)) #initial guess for the vector of decision variables
    u = forward(x_init)

#ubar is the location of the actual source
ubar = source_term(np.random.random_sample()*n, np.random.random_sample()*n, 2, sigma, W)
control = Control(u)


 #our functional with chosen source ubar (how can we had multiple ubars?
J = Functional((0.5*inner(u-ubar, u-ubar))*dx)


#bounds for each of the decision variables
lb = 0
ub = 1

rf = ReducedFunctional(J, control)

problem = MinimizationProblem(rf, bounds=(lb, ub))
parameters = {"acceptable_tol": 1.0e-3, "maximum_iterations": 100}
solver = IPOPTSolver(problem, parameters=parameters)
w_opt = solver.solve()
xdmf_filename = XDMFFile(mpi_comm_world(), "output/final_solution.xdmf")
xdmf_filename.write(w_opt)



enter code here