Switching a coefficient of term in the governing equations from 1 to 0 after some certain time

67 views
Skip to first unread message

Abishek Sarkar

unread,
Oct 25, 2024, 3:51:59 AMOct 25
to Dedalus Users
Dear all,

I am trying to switch on and switch off the forcing term in the governing equations (momentum equation) after a certain time (solver.sim_time > 0.5).
I initially define,
iForcing = 1
# Problem
problem = d3.IVP([u, p, tau_p], namespace=locals())
problem.add_equation("dt(u) + (1/rho_f)*grad(p) - nu*lap(u) = - u@grad(u) + iBody*fb + iForcing*Fw")
#problem.add_equation("dt(s) - D*lap(s) = - u@grad(s)")
problem.add_equation("div(u) + tau_p = 0")            # Used to put the boundary condition on pressure (Pressure Poisson's equation)
problem.add_equation("integ(p) = 0") # Pressure gauge

Say, I want to switch off the turbulent forcing in the momentum equation after t = 0.5 s, so I define a certain if-else condition,
# temporal variation of the swithces
if (solver. sim_time   < 0.5):
    iForcing = 1
else:
    iForcing = 0
But this condition does not seem to work. The governing equations are taking iForcing =1 , only. I even tried to define the if-else condition in the main loop, but it still doesn't seem to work.
Dear developers and fellow users, please help me out to tackle this issue of switching on and off the forcing term in the governing equations.

Thank you in advance.

Sinceely
Abishek Sarkar
 

Adrian Fraser

unread,
Oct 25, 2024, 10:15:59 AMOct 25
to Dedalus Users
Hi Abishek,

Have you taken a look at the documentation / posts on this list to do with GeneralFunctions? I'm thinking that's how you'd want to implement this feature.

Best,
Adrian

Curtis Saxton

unread,
Oct 28, 2024, 10:58:31 AMOct 28
to dedalu...@googlegroups.com
Hi Abishek,

I sometimes use a cruder trick to activate or deactivate a term at a certain time in the simulation. 

Run a simulation with one setting hardwired as constant, but stop at the simulation time when the change was due to occur (t=0.5 in your example). 

Edit the Dedalus code to change the relevant parameter or flag.  Hardwire the new value in your equations explicitly.  Restart the simulation using the last checkpoint file that was saved from the first stage simulation.

Would such a simple brute-force alteration be acceptable for your modelling?

Hope this helps,
Curtis Saxton.



--
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 visit https://groups.google.com/d/msgid/dedalus-users/350cd3c4-3084-4ff9-abb1-8230cd483118n%40googlegroups.com.

Miguel Beneitez

unread,
Oct 28, 2024, 10:58:31 AMOct 28
to dedalu...@googlegroups.com
Hi Abishek,

You can also do this directly by modifying the forcing (see below). This is probably not optimal but should work (I'm assuming here your forcing is constant and that it doesn't have a closed time dependent form). Also, I'd suggest that you maybe do this like: 
Fw['g'] = something # Initiliaze outside of the loop

# In loop
if (solver. sim_time  == 0.5): # so that it only goes in once
    Fw['g'] = 0.0

As this will save having to effectively assign the right hand side several times. If your forcing has a closed form or maybe shouldn't be updated like this on runtime, you can define an auxiliary field for iForcing, and then set it to be iForcing['g']=0/1 whenever you want. You carry one more field in memory but hopefully this is ok.

I hope this helps!
Miguel

See below for the suggested modification to your code.

# Problem
problem = d3.IVP([u, p, tau_p], namespace=locals())
problem.add_equation("dt(u) + (1/rho_f)*grad(p) - nu*lap(u) = - u@grad(u) + iBody*fb + Fw")

#problem.add_equation("dt(s) - D*lap(s) = - u@grad(s)")
problem.add_equation("div(u) + tau_p = 0")            # Used to put the boundary condition on pressure (Pressure Poisson's equation)
problem.add_equation("integ(p) = 0") # Pressure gauge

if (solver. sim_time   < 0.5):
    Fw['g'] = whatever you want it to be
else:
    Fw['g'] = whatever you want it to be

--

Abishek Sarkar

unread,
Oct 29, 2024, 10:01:47 AMOct 29
to Dedalus Users
Dear all,

Thank you so much for all your kind suggestions. I have looked upon the suggestion by Miguel and it gave me idea of putting a condition on the forcing term itself, rather than its coefficient, and it seems to work fine.

Thank you so much once again to all of you. And thank you so much to Miguel, as your suggestion helped me to solve the problem.

Sincerely
Abishek Sarkar

Reply all
Reply to author
Forward
0 new messages