External spatial and time dependent forcing

[Shane Kissinger]

Hello!

I'm still learning how to use Dedalus, so I apologize if there is an easy answer to this. I've searched on this group, but I wasn't able to find a post I could use. 

I'm trying to generate training data for the 1D burgers equation and want to add a randomized force term f(x,t). I have no issues adding a term f(t), but am struggling to implement a spatially and temporal dependent term. It's a snapshot of my overall code, but I've been using the below code:

################################ Bases #########################################

xcoord = d3.Coordinate('x')

dist = d3.Distributor(xcoord, dtype=dtype)

xbasis = d3.RealFourier(xcoord, size=grid, bounds=(x_min, x_max), dealias=dealias)

################################ Fields ########################################

u = dist.Field(name='u', bases=xbasis)

t = dist.Field(name='t')

############################# External Forcing #################################

seed = 0

k_min = 1

k_max = 3

rs = np.random.RandomState(seed)

a = 1 + rs.normal(0, 0.5)

omega = rs.uniform(-0.4, 0.4)

k_values = np.arange(k_min, k_max + 1)

k = rs.choice(k_values)

phi = rs.uniform(0, 2 * np.pi)

period = 2*np.pi

f = lambda x, t: a*np.sin(omega*t + (2 * np.pi * k * x) / period)

ext_force = ""

if forcing:

ext_force += " + f(x,t)" # Not sure what to put here either.

############################## Subsitutions ####################################

dx = lambda A: d3.Differentiate(A, xcoord)

################################ Problem #######################################

problem = d3.IVP([u], namespace=locals(), time='t')

equation = "dt(u) - eta*dx(dx(u)) = - u*dx(u)" + ext_force

problem.add_equation(equation)

How should I define the forcing function f(x,t)? It's just a sine function with two different parameters. 

 - Shane 

[Daniel Lecoanet]

Hi,

You want to use a GeneralFunction. You can search the users list for previous questions about how to implement a GeneralFunction.

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/5525de06-01c5-40a2-8e4e-d15ec884ba48n%40googlegroups.com.

[Shane Kissinger]

Hi Daniel,

I appreciate the response! I did see some prior posts about using a general function, and saw it was added to d3. However, I was confused when implementing it on what the “field” value should be if I am trying to access both x and t. 

Specifically, I’m talking about the implementation in the docs here

https://dedalus-project.readthedocs.io/en/latest/pages/general_functions.html

I think this might be stemming from a misunderstanding I have on how you can access the spatial coordinates to feed into the general function (as well as with t). The prior post I found was inputting random forcing via a pre-computed dataset. 

Sincerely,

Shane 

I’m 

[Shane Kissinger]

Hello!

I tried to implement a general function, though I'm not sure if this is correct. My goal is to have a randomized f(x,t) on the RHS that is a sin function dependent on t and on the x position within the 1d grid. The animations I'm running from time-stepping don't seem to be correct, though I could be mistaken. 

################################ Bases #########################################

xcoord = d3.Coordinate('x')

dist = d3.Distributor(xcoord, dtype=dtype)

xbasis = d3.RealFourier(xcoord, size=grid, bounds=(x_min, x_max), dealias=dealias)

domain = core.domain.Domain(dist=dist, bases=[xbasis])

################################ Fields ########################################

u = dist.Field(name='u', bases=xbasis)

t = dist.Field(name='t')

############################# External Forcing #################################

seed = 0

k_min = 1

k_max = 3

rs = np.random.RandomState(seed)

a = rs.normal(0, 0.5)

omega = rs.uniform(-0.4, 0.4)

k_values = np.arange(k_min, k_max + 1)

k = rs.choice(k_values)

phi = rs.uniform(0, 2 * np.pi)

period = 2*np.pi

def _force(tfield):

x = dist.local_grid(xbasis)

grid_dim = np.shape(x)

t = np.full(grid_dim, tfield.data)

return a*np.sin(omega*t + (2 * np.pi * k * x) / period)

def force(tfield):

return operators.GeneralFunction(

dist=dist,

domain=domain,

layout = 'g',

func = _force,

args = [tfield],

dtype=np.float64,

tensorsig=()

)

############################## Subsitutions ####################################

dx = lambda A: d3.Differentiate(A, xcoord)

################################ Problem #######################################

problem = d3.IVP([u], namespace=locals(), time='t')

equation = "dt(u) - eta*dx(dx(u)) = - u*dx(u) + force(t)"

problem.add_equation(equation)

############################ Initial Conditions ################################

x = dist.local_grid(xbasis)

u['g'] = a*np.sin((2 * np.pi * k * x) / period)

Sincerely, 

Shane