Field-Dependent Forcing Functions

[hitch...@gmail.com]

Hi,

I'm a student at Oxford and I've just started using Dedalus to study convection in mushy layers using the Rayleigh-Bernard example as a starting point. I'm having a bit of a problem setting up a forcing function that I hope someone might be able to help with though.

My liquid fraction is given by the function:
LiquidFraction = LiquidusGradient * BulkSalinity/( LiquidusGradient + InitialLiquidusTemperature - (BaseStateTemperature + PerturbationTemperature) )

I've been able to set this up as a standard equation but when the LiquidFraction goes above one, my equations break down and I get a numerical instability. To get around this, I want to cap this function so that it can't go above one.

(If there's a simple way to do this within the add_equation() format, stop me here!)

The solution I'm trying at the moment is to modify this solution:
https://groups.google.com/d/msg/dedalus-users/BqTjYZzqHHw/3iiohQHhh3sJ

to calculate the liquid fraction in a function defined elsewhere.

My code looks roughly like this:

# Liquid fraction function
def get_LqdFrac(solver, LqdusGrad, LqdusInf):
LqdFrac = LqdusGrad * solver.state['SalB']['g'] / (LqdusGrad + LqdusInf - (solver.state['Temp']['g'] + 1)*LqdusInf )
return LqdFrac

# Problem declaration
problem = de.IVP(domain, variables=['Temp','Temp_F','SalB','LqdF',...])
problem.meta['Temp','Temp_F','SalB',...]['z']['dirichlet'] = True

# Base-state temperature declaration
Temp_B = domain.new_field(name='Temp_B')
Temp_B['g'] = interp1d(Profiles['Z'],Profiles['Temp_B'],'cubic')(z)
Temp_B.meta['x']['constant'] = True
problem.parameters['Temp_B'] = Temp_B

# Forcing function declaration
LqdFrac_func = GeneralFunction(domain,'g',get_LqdFrac,args=[])
problem.parameters['LqdFrac_func'] = LqdFrac_func

# "Full" temperature equation and liquid fraction equation
problem.add_equation("Temp_F = Temp + Temp_B")
problem.add_equation("LqdF2 = LqdFrac_func")

# Adding function arguments
LqdFrac_func.args = [solver, Settings['LqdusGrad'], Settings['LqdusInf'] ]
LqdFrac_func.original_args = [solver, Settings['LqdusGrad'], Settings['LqdusInf'] ]

(I've missed out all the dynamical equations and other parts since they seem to be working fine but can include them if needed.)

This example above runs fine and produces sensible results for all the various fields, including Temp_F. My problem is that I need to use the full temperature (Temp_F) in my liquid fraction function rather than the perturbation temperature (Temp) but when I make this simple substitution I get the following error:

LqdFrac = LqdusGrad * solver.state['SalB']['g'] / (LqdusGrad + LqdusInf - (solver.state['Temp_F']['g'] + 1)*LqdusInf )
ValueError: operands could not be broadcast together with shapes (48,6) (32,4)

Something about swapping from "Temp" to "Temp_F" must be upsetting the code somehow but I don't understand what could be causing this problem when Temp_F seems to be calculated properly... I've tried substituting other fields in the space of Temp/Temp_F and they seem to run just fine, just like it does with Temp...

(Temp_B is also used in the dynamical equations and there are no issues there)

Any suggestions you guys have would be greatly appreciated!

Thanks,
Joe

[Daniel Lecoanet]

Hi Joe,

The problem is that that some of your fields are on the dealiased grid (size 48x6), whereas other fields are on the normal grid (32x4).  To fix this, before you calculate LqdFrac, you can add

solver.state['Temp_F'].set_scales(domain.dealias)

solver.state['SalB'].set_scales(domain.dealias)

Hope that helps,

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-users+unsubscribe@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/6a4fbadf-c38d-488c-b046-cdae96c6c8f5%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[hitch...@gmail.com]

> To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.

>
> To post to this group, send email to dedalu...@googlegroups.com.
>
> To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/6a4fbadf-c38d-488c-b046-cdae96c6c8f5%40googlegroups.com.
>
> For more options, visit https://groups.google.com/d/optout.

Thank you! That seems to have completely fixed the problem :)

[hitch...@gmail.com]