{"action":"created","issue":{"url":"https://api.github.com/repos/spectralDNS/shenfun/issues/18","repository_url":"https://api.github.com/repos/spectralDNS/shenfun","labels_url":"https://api.github.com/repos/spectralDNS/shenfun/issues/18/labels{/name}","comments_url":"https://api.github.com/repos/spectralDNS/shenfun/issues/18/comments","events_url":"https://api.github.com/repos/spectralDNS/shenfun/issues/18/events","html_url":"https://github.com/spectralDNS/shenfun/issues/18","id":357207561,"node_id":"MDU6SXNzdWUzNTcyMDc1NjE=","number":18,"title":"2D Rayleigh-Benard in vorticity-streamfunction formulation","user":{"login":"loiseaujc","id":1898398,"node_id":"MDQ6VXNlcjE4OTgzOTg=","avatar_url":"https://avatars1.githubusercontent.com/u/1898398?v=4","gravatar_id":"","url":"https://api.github.com/users/loiseaujc","html_url":"https://github.com/loiseaujc","followers_url":"https://api.github.com/users/loiseaujc/followers","following_url":"https://api.github.com/users/loiseaujc/following{/other_user}","gists_url":"https://api.github.com/users/loiseaujc/gists{/gist_id}","starred_url":"https://api.github.com/users/loiseaujc/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/loiseaujc/subscriptions","organizations_url":"https://api.github.com/users/loiseaujc/orgs","repos_url":"https://api.github.com/users/loiseaujc/repos","events_url":"https://api.github.com/users/loiseaujc/events{/privacy}","received_events_url":"https://api.github.com/users/loiseaujc/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":10,"created_at":"2018-09-05T12:33:50Z","updated_at":"2018-09-07T15:41:24Z","closed_at":null,"author_association":"NONE","body":"Hi Mikael.\r\n\r\nBased on shenfun, I have been able to implement a fairly efficient solver for two-dimensional decaying turbulence in a doubly-periodic domain which I eventually plan to use for part of my lectures. I am now willing to use shenfun again to implement a two-dimensional Rayleigh-Bénard solver using a vorticity-streamfunction formulation and I have a few questions. Before getting there, here are the equations I am interested in :\r\n\r\n<a href=\"https://www.codecogs.com/eqnedit.php?latex=\\begin{aligned}&space;\\frac{\\partial&space;\\omega}{\\partial&space;t}&space;&plus;&space;J(\\psi,&space;\\omega)&space;&&space;=&space;-Ra&space;Pr&space;\\frac{\\partial&space;\\theta}{\\partial&space;x}&space;&plus;&space;Pr&space;\\nabla^2&space;\\omega&space;\\\\&space;\\frac{\\partial&space;\\theta}{\\partial&space;t}&space;&plus;&space;J(\\psi,&space;\\theta)&space;&&space;=&space;\\frac{\\partial&space;\\psi}{\\partial&space;x}&space;&plus;&space;\\nabla^2&space;\\theta&space;\\\\&space;-\\nabla^2&space;\\psi&space;&&space;=&space;\\omega&space;\\end{aligned}\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?\\begin{aligned}&space;\\frac{\\partial&space;\\omega}{\\partial&space;t}&space;&plus;&space;J(\\psi,&space;\\omega)&space;&&space;=&space;-Ra&space;Pr&space;\\frac{\\partial&space;\\theta}{\\partial&space;x}&space;&plus;&space;Pr&space;\\nabla^2&space;\\omega&space;\\\\&space;\\frac{\\partial&space;\\theta}{\\partial&space;t}&space;&plus;&space;J(\\psi,&space;\\theta)&space;&&space;=&space;\\frac{\\partial&space;\\psi}{\\partial&space;x}&space;&plus;&space;\\nabla^2&space;\\theta&space;\\\\&space;-\\nabla^2&space;\\psi&space;&&space;=&space;\\omega&space;\\end{aligned}\" title=\"\\begin{aligned} \\frac{\\partial \\omega}{\\partial t} + J(\\psi, \\omega) & = -Ra Pr \\frac{\\partial \\theta}{\\partial x} + Pr \\nabla^2 \\omega \\\\ \\frac{\\partial \\theta}{\\partial t} + J(\\psi, \\theta) & = \\frac{\\partial \\psi}{\\partial x} + \\nabla^2 \\theta \\\\ -\\nabla^2 \\psi & = \\omega \\end{aligned}\" /></a>\r\n\r\nwhere <a href=\"https://www.codecogs.com/eqnedit.php?latex=\\omega\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?\\omega\" title=\"\\omega\" /></a> is the out-of-plane vorticity, <a href=\"https://www.codecogs.com/eqnedit.php?latex=\\psi\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?\\psi\" title=\"\\psi\" /></a> is the streamfunction and <a href=\"https://www.codecogs.com/eqnedit.php?latex=\\theta\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?\\theta\" title=\"\\theta\" /></a> is the deviation from the linear temperature profile. The set of boundary conditions needed to close the system are the following :\r\n- Homogeneous dirichlet conditions for <a href=\"https://www.codecogs.com/eqnedit.php?latex=\\theta\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?\\theta\" title=\"\\theta\" /></a>  on both walls\r\n- Bihamonic boundary conditions for  <a href=\"https://www.codecogs.com/eqnedit.php?latex=\\psi\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?\\psi\" title=\"\\psi\" /></a>  on both walls as well.\r\n\r\nFinally, all variables are periodic in the horizontal direction. The problem we consider is thus the canonical two-dimensional Rayleigh-Bénard convection with solid isothermal walls.\r\n\r\nAs far as my understanding of shenfun goes, here is basically what I need to do then :\r\n- I can use the same Fourier basis for all three variables, i.e.\r\n`F = Basis(N[1], family='F', dtype='D', domain=(-np.pi, np.pi))`\r\n\r\n- I then need to define the appropriate Chebyshev basis for each variable as\r\n`C_vorticity = Basis(N[0], family='C', domain=(0, 1))`\r\n`C_temperature = Basis(N[0], family='C', domain=(0, 1), bc='dirichlet')`\r\n`C_stream = Basis(N[0], family='C', domain=(0, 1), bc='biharmonic')`\r\n\r\n- Given these different bases, I can then define the appropriate TensorProductSpaces as\r\n`V_vorticity = TensorProductSpace(comm, (C_vorticity, F) **{'planner_effort': 'FFTW_PATIENT'}, axes=(0, 1))`\r\n`V_temperature = TensorProductSpace(comm, (C_temperature, F) **{'planner_effort': 'FFTW_PATIENT'}, axes=(0, 1))`\r\n`V_stream = TensorProductSpace(comm, (C_stream, F) **{'planner_effort': 'FFTW_PATIENT'}, axes=(0, 1))`\r\n\r\n- Finally, I can create the VectorProductSpace for my state vector <a href=\"https://www.codecogs.com/eqnedit.php?latex=Q&space;=&space;\\begin{bmatrix}&space;\\omega&space;&&space;\\theta&space;\\end{bmatrix}^T\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?Q&space;=&space;\\begin{bmatrix}&space;\\omega&space;&&space;\\theta&space;\\end{bmatrix}^T\" title=\"Q = \\begin{bmatrix} \\omega & \\theta \\end{bmatrix}^T\" /></a> as\r\n`Q = VectorTensorProductSpace([V_vorticity, V_temperature])`\r\n\r\nIf all of this is correct, I then have two questions, one related to the `LinearRHS` function, and the other to the `NonlinearRHS` ones. First, about `LinearRHS` : given my VectorTensorProductSpace Q, how can I create the linear operator\r\n\r\n<a href=\"https://www.codecogs.com/eqnedit.php?latex=\\mathcal{L}&space;=&space;\\begin{bmatrix}&space;Pr&space;\\nabla^2&space;&&space;0&space;\\\\&space;0&space;&&space;\\nabla^2&space;\\end{bmatrix}\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?\\mathcal{L}&space;=&space;\\begin{bmatrix}&space;Pr&space;\\nabla^2&space;&&space;0&space;\\\\&space;0&space;&&space;\\nabla^2&space;\\end{bmatrix}\" title=\"\\mathcal{L} = \\begin{bmatrix} Pr \\nabla^2 & 0 \\\\ 0 & \\nabla^2 \\end{bmatrix}\" /></a>\r\n\r\n( where I assume that I treat <a href=\"https://www.codecogs.com/eqnedit.php?latex=-Ra&space;Pr&space;\\frac{\\partial&space;\\theta}{\\partial&space;x}\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?-Ra&space;Pr&space;\\frac{\\partial&space;\\theta}{\\partial&space;x}\" title=\"-Ra Pr \\frac{\\partial \\theta}{\\partial x}\" /></a> and <a href=\"https://www.codecogs.com/eqnedit.php?latex=\\frac{\\partial&space;\\psi}{\\partial&space;x}\" target=\"_blank\"><img src=\"https://latex.codecogs.com/gif.latex?\\frac{\\partial&space;\\psi}{\\partial&space;x}\" title=\"\\frac{\\partial \\psi}{\\partial x}\" /></a> explicitly in the `NonlinearRHS` function)?\r\n\r\nIf creating such a block-diagonal linear operator is currently doable in shenfun, which type of linear solver does shenfun then uses to invert it considering the two different TensorProductSpaces used?\r\n\r\nThanks a lot,\r\nJC\r\n\r\n\r\nEDIT : I may have mixed some signs, but you get the points anyway."},"comment":{"url":"https://api.github.com/repos/spectralDNS/shenfun/issues/comments/419480190","html_url":"https://github.com/spectralDNS/shenfun/issues/18#issuecomment-419480190","issue_url":"https://api.github.com/repos/spectralDNS/shenfun/issues/18","id":419480190,"node_id":"MDEyOklzc3VlQ29tbWVudDQxOTQ4MDE5MA==","user":{"login":"loiseaujc","id":1898398,"node_id":"MDQ6VXNlcjE4OTgzOTg=","avatar_url":"https://avatars1.githubusercontent.com/u/1898398?v=4","gravatar_id":"","url":"https://api.github.com/users/loiseaujc","html_url":"https://github.com/loiseaujc","followers_url":"https://api.github.com/users/loiseaujc/followers","following_url":"https://api.github.com/users/loiseaujc/following{/other_user}","gists_url":"https://api.github.com/users/loiseaujc/gists{/gist_id}","starred_url":"https://api.github.com/users/loiseaujc/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/loiseaujc/subscriptions","organizations_url":"https://api.github.com/users/loiseaujc/orgs","repos_url":"https://api.github.com/users/loiseaujc/repos","events_url":"https://api.github.com/users/loiseaujc/events{/privacy}","received_events_url":"https://api.github.com/users/loiseaujc/received_events","type":"User","site_admin":false},"created_at":"2018-09-07T15:41:24Z","updated_at":"2018-09-07T15:41:24Z","author_association":"NONE","body":"I ran the same test with the 1D Poisson solver and homogeneous dirichlet boundary conditions. Here is the output:\r\n\r\n```\r\nIn [72]: run dirichlet_poisson1D.py 16 chebyshev\r\nError=5.2901900972253749e-06\r\n\r\nIn [73]: f_bis = f_hat.copy()\r\n\r\nIn [74]: f_bis.fill(0)\r\n\r\nIn [75]: f_bis = A.matvec(u_hat, f_bis)\r\n\r\nIn [76]: f_bis-f_hat\r\nOut[76]: \r\nFunction([ 0.00000000e+00,  0.00000000e+00, -1.57772181e-30,\r\n           0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\r\n           0.00000000e+00,  0.00000000e+00, -1.57772181e-30,\r\n           0.00000000e+00, -9.86076132e-32,  0.00000000e+00,\r\n           0.00000000e+00,  0.00000000e+00, -1.31703830e+01,\r\n           1.31703830e+01])\r\n\r\n```\r\nIt indeed seems like the problem comes from the boundary conditions, whether they are homogeneous or not."}}
{"actor":{"display_login":"loiseaujc"}}