Apply a filter to solution between iterations

1

Hi

I am trying to solve “Convection in the Earth’s Mantle” provided in the FEniCS book. Currently I don’t know how to handle the filter for the composition “filter_properties” developed by Lenardic and Kaula (1993)

I know you can access data from the solution of a PDE using phi.vector().get_local() as a numpy array, which is easy to manipulate, but I could not figure out how to convert the new values back to a function and use it for the next iteration.

The problem iterates over many timesteps, so the filtered \phi^k is necessary as a function to compute \phi^{k+1}

the filter bounds the value of the composition \phi between 0 and 1

What the filter is expected to do
(1) compute the sum of all values of composition \phi, S0
(2) Find Minimal value below 0 , assign to \phi_{min}
(3) Find Maximal value over 1, assign to \phi_{max}
(4) Assign to \phi \leq |\phi_{min}| value 0
(5) Assign to \phi \geq \phi_{max} value 1
(6) Compute the sum of all values of composition \phi, S1
(7) Compute the number num of 0 < \phi< 1
(8) Add dist = (S1 - S0)/num to all 0 < \phi < 1

Code example

  def filter_properties(phi):
        phi_numpy = phi.vector().get_local()
        S0 = np.sum(phi_numpy)
        phi_min = np.min(phi_numpy)
        phi_min = np.min([phi_min, 0])# in case phi_min > 0
        phi_max = np.max(phi)
        phi_max = np.max([phi_max, 1])# in case phi_max < 1

        for j in np.arange(len(phi_numpy)):
            if phi_numpy[j] < np.abs(phi_min):
                phi_numpy[j] = 0
            elif phi_numpy[j] > 2 - phi_max:
                phi_numpy[j] = 1

        S1 = np.sum(phi_numpy)
        num = 0

        for j in np.arange(len(phi_numpy)):
            if not(phi_numpy[j]==0 or phi_numpy[j] == 1):
                num += 1

        dist  = (S1 - S0)/num

        for j in np.arange(len(phi_numpy)):
            if not(phi[j]_numpy==0 or phi[j]_numpy == 1):
                phi_numpy[j] += dist

       #convert phi_numpy back to Function, How?
       return phi_Function

Is there any way to assign the filtered new values in \phi_{numpy} to reemplace the values in the function \phi
Thanks in advance

2

This paper describes the numerical method employed and provides example scripts for FEniCS in the supplementary material.

Sadly I cannot access the supplement at the moment since I don’t have access to the journal. But I remember it being written with FEniCS ~1.5. It will take a little work to update this to more recent versions such as 2019.2.0, but those differences will be mostly cosmetic.

1 Like
3

Thankyou very much, I could access the supplement. They have a script just for filters alone.
I will post the solution when I disentangle the code, for future coders that mayhavethe same question

1 Like
4

Here is an example

def Filter(phi):
values = phi.vector()[:] #recover values of phi as numpy array
# do the filtering algorithm to values
#--------------------------------------------

filtered = Function(phi.function_space) #  generate a new function
filtered.vector()[:] = values

return filtered
5

Hi, I am solving a similar problem, but it seems to me that the filter does not help with iterations. Did he help you?

6

Yes I managed to filter through each iteration. When you call the function in should be something like this

function.assign(filter(function)) # this inside your loop

2 Likes
7

I solve the problem of mixing two streams, my mesh is large enough and the Peclet number is high enough, but I noticed that discontinuous Galerkin methods are more numerically stable, but my concentration without a filter is greater than 1. Do you solve similar problems? Thank you for your answers.

8

Yes, I used the filter to bound composition (which can be thought of a dimensionless concentration) between 0 and 1.

I suggest you looking “Convection in Earth’s Mantle” in the FEniCS book, there you can find the refence for the filter used and it’s propertioes.

Also you should check your CFL condition if you are using an explicit method of discretization, that can be the source of non-physical concentrations.

1 Like
9

Yes, I follow the CFL, I try to choose the time step for the CFL = 0.25.
have you used such a filter? Can I borrow it for my solution?

10 
def filter_properties(phi):

    # filtro para acotar los valores del escalar phi entre 0 y 1
    # retorna funcion filtrada
    # extraer copia de phi como un numpy array, para poder manipular los datos

    vector = phi.vector()[:] 
    space = phi.function_space()

    # valores minimos y maximos

    vectorMin = np.min(vector)
    vectorMax = np.max(vector)

    # definimos la función filtrada, para ser retornada
    filtered = Function(space)
    S0 = np.sum(vector)
    n = len(vector)

    # filtrar valores bajo 0 a 0 y sobre 1 a 1

    vector = np.where(vector < 0, 0, vector)
    vector = np.where(vector > 1, 1, vector)

    # filtrar valoresmenos a abs(min) a 0 y sobre 2 - max a 1

    vector = np.where(vector < np.abs(vectorMin), 0, vector)
    vector = np.where(vector > (2 - vectorMax), 1, vector)

    # contamos los valores en el rango (0,1)

    num = 0
    for j in np.arange(n):

        if not(vector[j] == 0 or vector[j] == 1):

            num += 1

    # si no hay error que distribuirse retorna

    if num == 0:

        filtered.vector()[:] = vector

        return filtered

    S1 = np.sum(vector)

    dist = (S0-S1)/num

    for j in np.arange(n):

        if (vector[j] > 0 and vector[j] < 1):

            vector[j] += dist

    # volvemos al espacio de funciones asignando los valores filtrados al vector funcion 

    # filtrada

    filtered.vector()[:] = vector

    return filtered