Solving BVP together with IVP

[Andrei Igoshev]

Dear all,

I have a problem with Dedalus 3. I try to solve two coupled equations: one is initial value problem and second is boundary value problem. The IVP is solved just fine, but in the case of the BVP I find different solutions every time I try to solve it. Essentially the previous solution is modified if I give it as a guess solution.

The equation for boundary value problem is attached to this email. I have also attached the figures showing two consecutive tries to solve the problem. Essentially the solution looks very similar, but absolute numbers keep increasing. The problem is related to presence of d\delta \mu / dr in the right hand side of the equation. If I decrease this term I find exactly the same solution every time. It looks like I cannot move this term to the left-hand side because it is a derivative over one particular direction.

Essential code for BVP is attached to this email. The complete code which solves two coupled equations is available here: https://github.com/ignotur/ambipolar/blob/main/ambipolar1.py

Best regards,

Andrei Igoshev

[Andrei Igoshev]

I guess my question can be formulated also this way:

I have an BVP equation in form:

lap(x) - b * dot(grad(x), er) = f(y)

Is it possible to use radial derivative of x in the left-hand side (i.e.  b * dot(grad(x), er) ) using Dedalus?

If I try to use it now, I get the following error:

Traceback (most recent call last):

  File "/Users/amtai/Dropbox/research/work_with_rainer/ambipolar_diffusion/new_implementation/bvp_separate.py", line 132, in <module>

    solver = problem.build_solver()

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/problems.py", line 255, in build_solver

    return self.solver_class(self, *args, **kw)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/solvers.py", line 239, in __init__

    self.subproblems = subsystems.build_subproblems(self, self.subsystems, ['L'])

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/subsystems.py", line 74, in build_subproblems

    subproblem.build_matrices(matrices)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/subsystems.py", line 423, in build_matrices

    eqn_blocks = eqn[name].expression_matrices(subproblem=self, vars=vars, ncc_cutoff=solver.ncc_cutoff, max_ncc_terms=solver.max_ncc_terms)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/arithmetic.py", line 183, in expression_matrices

    arg_matrices = arg.expression_matrices(subproblem, vars, **kw)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/arithmetic.py", line 183, in expression_matrices

    arg_matrices = arg.expression_matrices(subproblem, vars, **kw)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/operators.py", line 711, in expression_matrices

    operator_mat = self.subproblem_matrix(subproblem)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/operators.py", line 2933, in subproblem_matrix

    factors[self.last_axis] = self.radial_matrix(regindex_in, regindex_out, ell)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/tools/cache.py", line 86, in __call__

    result = self.function(*args, **kw)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/operators.py", line 3744, in radial_matrix

    return self._radial_matrix(radial_basis, regtotal, ell)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/tools/cache.py", line 86, in __call__

    result = self.function(*args, **kw)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/operators.py", line 3751, in _radial_matrix

    return radial_basis.operator_matrix('L', ell, regtotal)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/tools/cache.py", line 86, in __call__

    result = self.function(*args, **kw)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/basis.py", line 3784, in operator_matrix

    size = self.n_size(l)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/tools/cache.py", line 86, in __call__

    result = self.function(*args, **kw)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/basis.py", line 3406, in n_size

    nmin = self._nmin(ell)

  File "/opt/anaconda3/envs/dedalus3/lib/python3.9/site-packages/dedalus/core/basis.py", line 3815, in _nmin

    return ell // 2

TypeError: unsupported operand type(s) for //: 'NoneType' and 'int'

[Evan Henry Anders]

Hi Andrei,

I'm not entirely sure what the solution is to your first problem. In response to your second question: yes, it is definitely possible to do "lap(x) - b * dot(grad(x), er) = f(y)" in dedalus. I'm attaching a script that does just that for the simplest case of f(y) = 0, and it runs successfully on my local laptop on commit c153f2ea6... on the master branch.

Perhaps try starting from here and see if you can build up to the full problem you're interested in? Also, can you take a look at my code here and see how it's different from what you've been trying and post that in a response (for people in the future who run into the same problem)?

Best,

Evan

--
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/0bc5f861-bf85-46e7-b5e4-93ebb9b7b569n%40googlegroups.com.

[Andrei Igoshev]

Dear Evan,

Thank you! Your example works for me. There is just one problem. 

I have a complicated radial profile as my b coefficient.

lap(x) - b (r) * dot(grad(x), er) = f(y)

Is it possible to use b(r) as a factor in the left-hand side of the equation? 

If I modify your script and add this radial profile I have got the same error:

python3 example_dedalus.py

Traceback (most recent call last):

  File "/Users/amtai/Dropbox/research/work_with_rainer/ambipolar_diffusion/new_implementation/example_dedalus.py", line 43, in <module>

You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/Ac0YMPv7Nh8/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CALW8q5Ob4Ztz%3DjfsjPEt-uq6s_PZLK%3DCZBDgTJ8DrCkvUok2Uw%40mail.gmail.com.

[Evan Henry Anders]

Hi Andrei,

Great! Yes, you can do non-constant coefficients like you've mentioned in d3 (e.g., using b(r) ), you just need to define the field for the nonconstant coefficient on the *radial* basis instead of on the full basis. I've slightly modified the code you sent to do this and it runs with b(r) = a*r^2 for me.

One thing to be careful of: b(r) needs to not violate the regularity of the Ball's radial basis. I think this mostly means that its derivative has to vanish at r = 0 (so dr(b)(r=0) must be 0), but there may be other things I'm not aware of.

Best,

Evan

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAMa8c%2BSaE7gc8%2B6Ki8UAXZV6O47bNNjgh3nTKZ8Qk3a_WP1n%2Bg%40mail.gmail.com.

[Andrei Igoshev]

Dear Evan,

Thank you very much, it solves my problem!

Best regards,

Andrei Igoshev