[![binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gl/rubendibattista%2Fjosiepy/master/?filepath=examples%2F002_Advection.ipynb) # Implementing a solver for the advection equation We're going to write all the necessary objects to solve the following equation: ```math \partial_t \varphi + \mathbf{u} \cdot \nabla \varphi = 0 ``` That, compared to our generalized model: ```math \partial_t \mathbf{q} + \nabla \cdot \left(\underline{\underline{\mathbf{F}}}(\mathbf{q}) + \underline{\underline{\underline{\underline{D(\mathbf{q})}}}} \cdot \nabla \mathbf{q} \right) + \underline{\underline{\underline{\mathbf{B}(\mathbf{q})}}} \cdot \nabla\mathbf{q} = \mathbf{s}(\mathbf{q}) ``` Gives `$ \mathbf{q} = \left\{\varphi\right\} $` `$ \mathbf{F} = \underline{\underline{\mathbf{0}}} $`, `$ \mathbf{D} = \underline{\underline{\underline{\underline{\mathbf{0}}}}} $` So our problem needs just to provide `$ \mathbf{B} $`. ```math B_{pqr} = u_r ``` ## `State` ```python from josie.state import State from josie.fields import Fields class AdvectionFields(Fields): phi = 0 class AdvectionState(State): fields = AdvectionFields ``` ## `Problem` We just need to implement the `$ \mathbf{B}(\mathbf{q}) $` operator ```python from josie.problem import Problem Problem?? ``` ```python from josie.mesh.cellset import CellSet, MeshCellSet class AdvectionProblem(Problem): V = np.array([1, 0]) DIMENSIONALITY = 2 ## 2D def B(self, cells: CellSet): values: AdvectionState = cells.values nx, ny, num_fields = values.shape fields = values.fields B = np.zeros((nx, ny, num_fields, num_fields, self.DIMENSIONALITY)) B[..., fields.phi, fields.phi, : ] = self.V return B ``` ## `Scheme` ![cell_scheme](data:image/svg+xml;base64,PD94bWwgdmVyc2lvbj0iMS4wIiBlbmNvZGluZz0iVVRGLTgiIHN0YW5kYWxvbmU9Im5vIj8+CjxzdmcKICAgeG1sbnM6ZGM9Imh0dHA6Ly9wdXJsLm9yZy9kYy9lbGVtZW50cy8xLjEvIgogICB4bWxuczpjYz0iaHR0cDovL2NyZWF0aXZlY29tbW9ucy5vcmcvbnMjIgogICB4bWxuczpyZGY9Imh0dHA6Ly93d3cudzMub3JnLzE5OTkvMDIvMjItcmRmLXN5bnRheC1ucyMiCiAgIHhtbG5zOnN2Zz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciCiAgIHhtbG5zPSJodHRwOi8vd3d3LnczLm9yZy8yMDAwL3N2ZyIKICAgeG1sbnM6c29kaXBvZGk9Imh0dHA6Ly9zb2RpcG9kaS5zb3VyY2Vmb3JnZS5uZXQvRFREL3NvZGlwb2RpLTAuZHRkIgogICB4bWxuczppbmtzY2FwZT0iaHR0cDovL3d3dy5pbmtzY2FwZS5vcmcvbmFtZXNwYWNlcy9pbmtzY2FwZSIKICAgaW5rc2NhcGU6dmVyc2lvbj0iMS4wcmMxICgwOTk2MGQ2LCAyMDIwLTA0LTA5KSIKICAgc29kaXBvZGk6ZG9jbmFtZT0iY2VsbF9zY2hlbWUuc3ZnIgogICB2aWV3Qm94PSIwIDAgMTk3LjMzMzMzIDE0NCIKICAgaGVpZ2h0PSIxNDQiCiAgIHdpZHRoPSIxOTcuMzMzMzMiCiAgIHhtbDpzcGFjZT0icHJlc2VydmUiCiAgIGlkPSJzdmcxMDQwIgogICB2ZXJzaW9uPSIxLjEiPjxtZXRhZGF0YQogICAgIGlkPSJtZXRhZGF0YTEwNDYiPjxyZGY6UkRGPjxjYzpXb3JrCiAgICAgICAgIHJkZjphYm91dD0iIj48ZGM6Zm9ybWF0PmltYWdlL3N2Zyt4bWw8L2RjOmZvcm1hdD48ZGM6dHlwZQogICAgICAgICAgIHJkZjpyZXNvdXJjZT0iaHR0cDovL3B1cmwub3JnL2RjL2RjbWl0eXBlL1N0aWxsSW1hZ2UiIC8+PGRjOnRpdGxlPjwvZGM6dGl0bGU+PC9jYzpXb3JrPjwvcmRmOlJERj48L21ldGFkYXRhPjxkZWZzCiAgICAgaWQ9ImRlZnMxMDQ0IiAvPjxzb2RpcG9kaTpuYW1lZHZpZXcKICAgICBpbmtzY2FwZTpjdXJyZW50LWxheWVyPSJnMTA0OCIKICAgICBpbmtzY2FwZTp3aW5kb3ctbWF4aW1pemVkPSIwIgogICAgIGlua3NjYXBlOndpbmRvdy15PSIyMyIKICAgICBpbmtzY2FwZTp3aW5kb3cteD0iMCIKICAgICBpbmtzY2FwZTpjeT0iNzIiCiAgICAgaW5rc2NhcGU6Y3g9Ijk4LjY2NjY2NCIKICAgICBpbmtzY2FwZTp6b29tPSIzLjYwODEwODIiCiAgICAgc2hvd2dyaWQ9ImZhbHNlIgogICAgIGlkPSJuYW1lZHZpZXcxMDQyIgogICAgIGlua3NjYXBlOndpbmRvdy1oZWlnaHQ9Ijk1MyIKICAgICBpbmtzY2FwZTp3aW5kb3ctd2lkdGg9IjE2ODAiCiAgICAgaW5rc2NhcGU6cGFnZXNoYWRvdz0iMiIKICAgICBpbmtzY2FwZTpwYWdlb3BhY2l0eT0iMCIKICAgICBndWlkZXRvbGVyYW5jZT0iMTAiCiAgICAgZ3JpZHRvbGVyYW5jZT0iMTAiCiAgICAgb2JqZWN0dG9sZXJhbmNlPSIxMCIKICAgICBib3JkZXJvcGFjaXR5PSIxIgogICAgIGJvcmRlcmNvbG9yPSIjNjY2NjY2IgogICAgIHBhZ2Vjb2xvcj0iI2ZmZmZmZiIgLz48ZwogICAgIHRyYW5zZm9ybT0ibWF0cml4KDEuMzMzMzMzMywwLDAsLTEuMzMzMzMzMywwLDE0NCkiCiAgICAgaW5rc2NhcGU6bGFiZWw9Imlua19leHRfWFhYWFhYIgogICAgIGlua3NjYXBlOmdyb3VwbW9kZT0ibGF5ZXIiCiAgICAgaWQ9ImcxMDQ4Ij48ZwogICAgICAgdHJhbnNmb3JtPSJzY2FsZSgwLjEpIgogICAgICAgaWQ9ImcxMDUwIj48cGF0aAogICAgICAgICBpbmtzY2FwZTpjb25uZWN0b3ItY3VydmF0dXJlPSIwIgogICAgICAgICBpZD0icGF0aDEwNTIiCiAgICAgICAgIHN0eWxlPSJmaWxsOm5vbmU7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjg7c3Ryb2tlLWxpbmVjYXA6YnV0dDtzdHJva2UtbGluZWpvaW46cm91bmQ7c3Ryb2tlLW1pdGVybGltaXQ6MTA7c3Ryb2tlLWRhc2hhcnJheTpub25lO3N0cm9rZS1vcGFjaXR5OjEiCiAgICAgICAgIGQ9Im0gODAzLjk4LDY0NS43MzggaCA2NDAgViA1LjczODI4IGwgLTQ4MCwxNTkuOTk5NzIgLTE2MCw0ODAiIC8+PHBhdGgKICAgICAgICAgaW5rc2NhcGU6Y29ubmVjdG9yLWN1cnZhdHVyZT0iMCIKICAgICAgICAgaWQ9InBhdGgxMDU0IgogICAgICAgICBzdHlsZT0iZmlsbDpub25lO3N0cm9rZTojMDAwMDAwO3N0cm9rZS13aWR0aDo4O3N0cm9rZS1saW5lY2FwOmJ1dHQ7c3Ryb2tlLWxpbmVqb2luOnJvdW5kO3N0cm9rZS1taXRlcmxpbWl0OjEwO3N0cm9rZS1kYXNoYXJyYXk6bm9uZTtzdHJva2Utb3BhY2l0eToxIgogICAgICAgICBkPSJNIDgwMy45OCw2NDQuMDIgSCAzLjk4MDQ3IEwgMTYzLjk4LDQuMDE5NTMgOTYzLjk4LDE2NC4wMiIgLz48cGF0aAogICAgICAgICBpbmtzY2FwZTpjb25uZWN0b3ItY3VydmF0dXJlPSIwIgogICAgICAgICBpZD0icGF0aDEwNTYiCiAgICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7ZmlsbC1vcGFjaXR5OjE7ZmlsbC1ydWxlOmV2ZW5vZGQ7c3Ryb2tlOm5vbmUiCiAgICAgICAgIGQ9Im0gNDkwLjU1MSwzNzIuNzcgYyAwLDI0LjAzMSAtMzYuMDIsMjQuMDMxIC0zNi4wMiwwIDAsLTIzLjk4MSAzNi4wMiwtMjMuOTgxIDM2LjAyLDAgeiIgLz48cGF0aAogICAgICAgICBpbmtzY2FwZTpjb25uZWN0b3ItY3VydmF0dXJlPSIwIgogICAgICAgICBpZD0icGF0aDEwNTgiCiAgICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7ZmlsbC1vcGFjaXR5OjE7ZmlsbC1ydWxlOmV2ZW5vZGQ7c3Ryb2tlOm5vbmUiCiAgICAgICAgIGQ9Im0gMTIwOS4xNCw0MDQuNDQ5IGMgMCwyMy45OTIgLTM2LjAyLDIzLjk5MiAtMzYuMDIsMCAwLC0yNC4wMTkgMzYuMDIsLTI0LjAxOSAzNi4wMiwwIHoiIC8+PHBhdGgKICAgICAgICAgaW5rc2NhcGU6Y29ubmVjdG9yLWN1cnZhdHVyZT0iMCIKICAgICAgICAgaWQ9InBhdGgxMDYwIgogICAgICAgICBzdHlsZT0iZmlsbDpub25lO3N0cm9rZTojMDAwMDAwO3N0cm9rZS13aWR0aDo0O3N0cm9rZS1saW5lY2FwOmJ1dHQ7c3Ryb2tlLWxpbmVqb2luOnJvdW5kO3N0cm9rZS1taXRlcmxpbWl0OjEwO3N0cm9rZS1kYXNoYXJyYXk6bm9uZTtzdHJva2Utb3BhY2l0eToxIgogICAgICAgICBkPSJNIDY4MC4xOTksMzM5LjY0OCAxMDg5LjQ5LDQ5Mi4xOTEiIC8+PHBhdGgKICAgICAgICAgaW5rc2NhcGU6Y29ubmVjdG9yLWN1cnZhdHVyZT0iMCIKICAgICAgICAgaWQ9InBhdGgxMDYyIgogICAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO2ZpbGwtb3BhY2l0eToxO2ZpbGwtcnVsZTpldmVub2RkO3N0cm9rZTpub25lIgogICAgICAgICBkPSJtIDEwODkuNDksNDkyLjE5MSAtNzMuNzUsLTIuNjIxIDE2LjI5LC00My42NzIgeiIgLz48cGF0aAogICAgICAgICBpbmtzY2FwZTpjb25uZWN0b3ItY3VydmF0dXJlPSIwIgogICAgICAgICBpZD0icGF0aDEwNjQiCiAgICAgICAgIHN0eWxlPSJmaWxsOm5vbmU7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjQ7c3Ryb2tlLWxpbmVjYXA6YnV0dDtzdHJva2UtbGluZWpvaW46cm91bmQ7c3Ryb2tlLW1pdGVybGltaXQ6MTA7c3Ryb2tlLWRhc2hhcnJheTpub25lO3N0cm9rZS1vcGFjaXR5OjEiCiAgICAgICAgIGQ9Im0gMTA4OS40OSw0OTIuMTkxIC03My43NSwtMi42MjEgMTYuMjksLTQzLjY3MiB6IiAvPjxnCiAgICAgICAgIHRyYW5zZm9ybT0ic2NhbGUoMTApIgogICAgICAgICBpZD0iZzEwNjYiPjx0ZXh0CiAgICAgICAgICAgaWQ9InRleHQxMDcwIgogICAgICAgICAgIHN0eWxlPSJmb250LXZhcmlhbnQ6bm9ybWFsO2ZvbnQtc2l6ZTo5Ljk2MjZweDtmb250LWZhbWlseTpmLTA7LWlua3NjYXBlLWZvbnQtc3BlY2lmaWNhdGlvbjpmLTAtMDt3cml0aW5nLW1vZGU6bHItdGI7ZmlsbDojMDAwMDAwO2ZpbGwtb3BhY2l0eToxO2ZpbGwtcnVsZTpub256ZXJvO3N0cm9rZTpub25lIgogICAgICAgICAgIHRyYW5zZm9ybT0ibWF0cml4KDEsMCwwLC0xLDkyLjk5ODQsNDguNjg1MikiPjx0c3BhbgogICAgICAgICAgICAgaWQ9InRzcGFuMTA2OCIKICAgICAgICAgICAgIHk9IjAiCiAgICAgICAgICAgICB4PSIwIj7LhjwvdHNwYW4+PC90ZXh0Pjx0ZXh0CiAgICAgICAgICAgaWQ9InRleHQxMDc0IgogICAgICAgICAgIHN0eWxlPSJmb250LXZhcmlhbnQ6bm9ybWFsO2ZvbnQtc2l6ZTo5Ljk2MjZweDtmb250LWZhbWlseTpmLTE7LWlua3NjYXBlLWZvbnQtc3BlY2lmaWNhdGlvbjpmLTEtMDt3cml0aW5nLW1vZGU6bHItdGI7ZmlsbDojMDAwMDAwO2ZpbGwtb3BhY2l0eToxO2ZpbGwtcnVsZTpub256ZXJvO3N0cm9rZTpub25lIgogICAgICAgICAgIHRyYW5zZm9ybT0ibWF0cml4KDEsMCwwLC0xLDkyLjMwNjIxMiw0OC41NDY5MTkpIj48dHNwYW4KICAgICAgICAgICAgIGlkPSJ0c3BhbjEwNzIiCiAgICAgICAgICAgICB5PSIwIgogICAgICAgICAgICAgeD0iMCI+bjwvdHNwYW4+PC90ZXh0PjwvZz48cGF0aAogICAgICAgICBpbmtzY2FwZTpjb25uZWN0b3ItY3VydmF0dXJlPSIwIgogICAgICAgICBpZD0icGF0aDEwNzYiCiAgICAgICAgIHN0eWxlPSJmaWxsOm5vbmU7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjIwO3N0cm9rZS1saW5lY2FwOmJ1dHQ7c3Ryb2tlLWxpbmVqb2luOnJvdW5kO3N0cm9rZS1taXRlcmxpbWl0OjEwO3N0cm9rZS1kYXNoYXJyYXk6bm9uZTtzdHJva2Utb3BhY2l0eToxIgogICAgICAgICBkPSJtIDgwMy45OCw2NDQuMDIgMTYwLC00ODAiIC8+PGcKICAgICAgICAgdHJhbnNmb3JtPSJzY2FsZSgxMCkiCiAgICAgICAgIGlkPSJnMTA3OCI+PHRleHQKICAgICAgICAgICBpZD0idGV4dDEwODIiCiAgICAgICAgICAgc3R5bGU9ImZvbnQtdmFyaWFudDpub3JtYWw7Zm9udC1zaXplOjkuOTYyNnB4O2ZvbnQtZmFtaWx5OmYtMjstaW5rc2NhcGUtZm9udC1zcGVjaWZpY2F0aW9uOmYtMi0wO3dyaXRpbmctbW9kZTpsci10YjtmaWxsOiMwMDAwMDA7ZmlsbC1vcGFjaXR5OjE7ZmlsbC1ydWxlOm5vbnplcm87c3Ryb2tlOm5vbmUiCiAgICAgICAgICAgdHJhbnNmb3JtPSJtYXRyaXgoMSwwLDAsLTEsOTUuNjY0NSwzMC42ODc5KSI+PHRzcGFuCiAgICAgICAgICAgICBpZD0idHNwYW4xMDgwIgogICAgICAgICAgICAgeT0iMCIKICAgICAgICAgICAgIHg9IjAiPlM8L3RzcGFuPjwvdGV4dD48dGV4dAogICAgICAgICAgIGlkPSJ0ZXh0MTA4NiIKICAgICAgICAgICBzdHlsZT0iZm9udC12YXJpYW50Om5vcm1hbDtmb250LXNpemU6Ni45NzM4cHg7Zm9udC1mYW1pbHk6Zi0zOy1pbmtzY2FwZS1mb250LXNwZWNpZmljYXRpb246Zi0zLTA7d3JpdGluZy1tb2RlOmxyLXRiO2ZpbGw6IzAwMDAwMDtmaWxsLW9wYWNpdHk6MTtmaWxsLXJ1bGU6bm9uemVybztzdHJva2U6bm9uZSIKICAgICAgICAgICB0cmFuc2Zvcm09Im1hdHJpeCgxLDAsMCwtMSwxMDEuNzczNDgsMjkuMTkzNzYpIj48dHNwYW4KICAgICAgICAgICAgIGlkPSJ0c3BhbjEwODQiCiAgICAgICAgICAgICB5PSIwIgogICAgICAgICAgICAgeD0iMCI+ZjwvdHNwYW4+PC90ZXh0PjwvZz48cGF0aAogICAgICAgICBpbmtzY2FwZTpjb25uZWN0b3ItY3VydmF0dXJlPSIwIgogICAgICAgICBpZD0icGF0aDEwODgiCiAgICAgICAgIHN0eWxlPSJmaWxsOm5vbmU7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjQ7c3Ryb2tlLWxpbmVjYXA6YnV0dDtzdHJva2UtbGluZWpvaW46cm91bmQ7c3Ryb2tlLW1pdGVybGltaXQ6MTA7c3Ryb2tlLWRhc2hhcnJheTpub25lO3N0cm9rZS1vcGFjaXR5OjEiCiAgICAgICAgIGQ9Ik0gNDc0LjIxOSwzNzIuNzcgMzY1LjM1Miw4OTAuODk4IiAvPjxwYXRoCiAgICAgICAgIGlua3NjYXBlOmNvbm5lY3Rvci1jdXJ2YXR1cmU9IjAiCiAgICAgICAgIGlkPSJwYXRoMTA5MCIKICAgICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtmaWxsLW9wYWNpdHk6MTtmaWxsLXJ1bGU6ZXZlbm9kZDtzdHJva2U6bm9uZSIKICAgICAgICAgZD0ibSA0NzQuMjE5LDM3Mi43NyA4LjQwMiw3My4zMiAtNDUuNjMzLC05LjYxIHoiIC8+PHBhdGgKICAgICAgICAgaW5rc2NhcGU6Y29ubmVjdG9yLWN1cnZhdHVyZT0iMCIKICAgICAgICAgaWQ9InBhdGgxMDkyIgogICAgICAgICBzdHlsZT0iZmlsbDpub25lO3N0cm9rZTojMDAwMDAwO3N0cm9rZS13aWR0aDo0O3N0cm9rZS1saW5lY2FwOmJ1dHQ7c3Ryb2tlLWxpbmVqb2luOnJvdW5kO3N0cm9rZS1taXRlcmxpbWl0OjEwO3N0cm9rZS1kYXNoYXJyYXk6bm9uZTtzdHJva2Utb3BhY2l0eToxIgogICAgICAgICBkPSJtIDQ3NC4yMTksMzcyLjc3IDguNDAyLDczLjMyIC00NS42MzMsLTkuNjEgeiIgLz48ZwogICAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEwKSIKICAgICAgICAgaWQ9ImcxMDk0Ij48dGV4dAogICAgICAgICAgIGlkPSJ0ZXh0MTA5OCIKICAgICAgICAgICBzdHlsZT0iZm9udC12YXJpYW50Om5vcm1hbDtmb250LXNpemU6OS45NjI2cHg7Zm9udC1mYW1pbHk6Zi00Oy1pbmtzY2FwZS1mb250LXNwZWNpZmljYXRpb246Zi00LTA7d3JpdGluZy1tb2RlOmxyLXRiO2ZpbGw6IzAwMDAwMDtmaWxsLW9wYWNpdHk6MTtmaWxsLXJ1bGU6bm9uemVybztzdHJva2U6bm9uZSIKICAgICAgICAgICB0cmFuc2Zvcm09Im1hdHJpeCgxLDAsMCwtMSwyNC4xODc5LDkwLjA3NDYpIj48dHNwYW4KICAgICAgICAgICAgIGlkPSJ0c3BhbjEwOTYiCiAgICAgICAgICAgICBzb2RpcG9kaTpyb2xlPSJsaW5lIgogICAgICAgICAgICAgeT0iMCIKICAgICAgICAgICAgIHg9IjAgNS4yMzAzNjQ4IDEwLjQ2MDczIDE1LjY5MTA5NSAyMC45MjE0NTkiPnZhbHVlPC90c3Bhbj48L3RleHQ+PC9nPjxwYXRoCiAgICAgICAgIGlua3NjYXBlOmNvbm5lY3Rvci1jdXJ2YXR1cmU9IjAiCiAgICAgICAgIGlkPSJwYXRoMTEwMCIKICAgICAgICAgc3R5bGU9ImZpbGw6bm9uZTtzdHJva2U6IzAwMDAwMDtzdHJva2Utd2lkdGg6NDtzdHJva2UtbGluZWNhcDpidXR0O3N0cm9rZS1saW5lam9pbjpyb3VuZDtzdHJva2UtbWl0ZXJsaW1pdDoxMDtzdHJva2UtZGFzaGFycmF5Om5vbmU7c3Ryb2tlLW9wYWNpdHk6MSIKICAgICAgICAgZD0ibSAxMTkxLjEzLDQwNC40NDkgLTMyLjY1LDU0MC42NzIiIC8+PHBhdGgKICAgICAgICAgaW5rc2NhcGU6Y29ubmVjdG9yLWN1cnZhdHVyZT0iMCIKICAgICAgICAgaWQ9InBhdGgxMTAyIgogICAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO2ZpbGwtb3BhY2l0eToxO2ZpbGwtcnVsZTpldmVub2RkO3N0cm9rZTpub25lIgogICAgICAgICBkPSJtIDExOTEuMTMsNDA0LjQ0OSAxOS4wNyw3MS4yODkgLTQ2LjU3LC0yLjgwOCB6IiAvPjxwYXRoCiAgICAgICAgIGlua3NjYXBlOmNvbm5lY3Rvci1jdXJ2YXR1cmU9IjAiCiAgICAgICAgIGlkPSJwYXRoMTEwNCIKICAgICAgICAgc3R5bGU9ImZpbGw6bm9uZTtzdHJva2U6IzAwMDAwMDtzdHJva2Utd2lkdGg6NDtzdHJva2UtbGluZWNhcDpidXR0O3N0cm9rZS1saW5lam9pbjpyb3VuZDtzdHJva2UtbWl0ZXJsaW1pdDoxMDtzdHJva2UtZGFzaGFycmF5Om5vbmU7c3Ryb2tlLW9wYWNpdHk6MSIKICAgICAgICAgZD0ibSAxMTkxLjEzLDQwNC40NDkgMTkuMDcsNzEuMjg5IC00Ni41NywtMi44MDggeiIgLz48ZwogICAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEwKSIKICAgICAgICAgaWQ9ImcxMTA2Ij48dGV4dAogICAgICAgICAgIGlkPSJ0ZXh0MTExMCIKICAgICAgICAgICBzdHlsZT0iZm9udC12YXJpYW50Om5vcm1hbDtmb250LXNpemU6OS45NjI2cHg7Zm9udC1mYW1pbHk6Zi00Oy1pbmtzY2FwZS1mb250LXNwZWNpZmljYXRpb246Zi00LTA7d3JpdGluZy1tb2RlOmxyLXRiO2ZpbGw6IzAwMDAwMDtmaWxsLW9wYWNpdHk6MTtmaWxsLXJ1bGU6bm9uemVybztzdHJva2U6bm9uZSIKICAgICAgICAgICB0cmFuc2Zvcm09Im1hdHJpeCgxLDAsMCwtMSw5Mi4xNzIzLDEwMC45MTIpIj48dHNwYW4KICAgICAgICAgICAgIGlkPSJ0c3BhbjExMDgiCiAgICAgICAgICAgICBzb2RpcG9kaTpyb2xlPSJsaW5lIgogICAgICAgICAgICAgeT0iMCIKICAgICAgICAgICAgIHg9IjAgNS4yMzAzNjQ4IDEwLjQ2MDczIDE1LjY5MTA5NSAyMC45MjE0NTkiPm5laWdoPC90c3Bhbj48L3RleHQ+PC9nPjxwYXRoCiAgICAgICAgIGlua3NjYXBlOmNvbm5lY3Rvci1jdXJ2YXR1cmU9IjAiCiAgICAgICAgIGlkPSJwYXRoMTExMiIKICAgICAgICAgc3R5bGU9ImZpbGw6bm9uZTtzdHJva2U6IzAwMDAwMDtzdHJva2Utd2lkdGg6My45ODtzdHJva2UtbGluZWNhcDpidXR0O3N0cm9rZS1saW5lam9pbjptaXRlcjtzdHJva2UtbWl0ZXJsaW1pdDoxMDtzdHJva2UtZGFzaGFycmF5Om5vbmU7c3Ryb2tlLW9wYWNpdHk6MSIKICAgICAgICAgZD0ibSAxMTg5LjUzLDEwMTEuMDkgaCAzMS4zNyIgLz48ZwogICAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEwKSIKICAgICAgICAgaWQ9ImcxMTE0Ij48dGV4dAogICAgICAgICAgIGlkPSJ0ZXh0MTExOCIKICAgICAgICAgICBzdHlsZT0iZm9udC12YXJpYW50Om5vcm1hbDtmb250LXNpemU6OS45NjI2cHg7Zm9udC1mYW1pbHk6Zi00Oy1pbmtzY2FwZS1mb250LXNwZWNpZmljYXRpb246Zi00LTA7d3JpdGluZy1tb2RlOmxyLXRiO2ZpbGw6IzAwMDAwMDtmaWxsLW9wYWNpdHk6MTtmaWxsLXJ1bGU6bm9uemVybztzdHJva2U6bm9uZSIKICAgICAgICAgICB0cmFuc2Zvcm09Im1hdHJpeCgxLDAsMCwtMSwxMjIuMDkxLDEwMC45MTIpIj48dHNwYW4KICAgICAgICAgICAgIGlkPSJ0c3BhbjExMTYiCiAgICAgICAgICAgICBzb2RpcG9kaTpyb2xlPSJsaW5lIgogICAgICAgICAgICAgeT0iMCIKICAgICAgICAgICAgIHg9IjAgNS4yMzAzNjQ4IDEwLjQ2MDczIDE1LjY5MTA5NSAyMC45MjE0NTkiPnZhbHVlPC90c3Bhbj48L3RleHQ+PC9nPjwvZz48L2c+PC9zdmc+Cg==) We just need to implement the ```math \mathbf{G}_{\frac{1}{2}}(\mathbf{q}) = \left| \mathbf{q}\hat{\mathbf{n}} \right|_f S_f \Rightarrow \left| \varphi \hat{\mathbf{n}}\right|_f S_f ``` ```python from josie.scheme.nonconservative import NonConservativeScheme NonConservativeScheme?? ``` ```python class Upwind(NonConservativeScheme): def G(self, cells: MeshCellSet, neighs: CellSet): V = self.problem.V # Get the normal velocities (that are equal everywhere) per each cell Vn = np.einsum("k,...kl->...l", V, neighs.normals[..., np.newaxis]) values_face = np.zeros_like(cells.values) # We use the cell value where Vn > 0 np.copyto(values_face, cells.values, where=(Vn>0)) # We use the neighbour value otherwise np.copyto(values_face, neighs.values, where=(Vn<=0)) # Multiply by the normal valuesn_face = np.einsum("...i,...j->...ij", values_face, neighs.normals) # Multiply by the surface G = valuesn_face * neighs.surfaces[..., np.newaxis, np.newaxis] return G ``` ## `Mesh` We generated a 1D mesh ```python from josie.boundary import Line left = Line([0, 0], [0, 1]) bottom = Line([0, 0], [1, 0]) right = Line([1, 0], [1, 1]) top = Line([0, 1], [1, 1]) ``` We apply periodic boundary condition along the x-axis (no BC on y-axis since it's a 1D simulation) ```python from josie.bc import make_periodic, Direction left, right = make_periodic(left, right, Direction.X) top.bc = None bottom.bc = None print(left, right, top, bottom) ``` ```python from josie.mesh import Mesh from josie.mesh.cell import SimpleCell mesh = Mesh(left, bottom, right, top, SimpleCell) mesh.interpolate(300, 1) mesh.generate() ``` ## `Solver` Let's assemble our scheme, that still needs the time scheme, and a CFL method ```python from josie.general.schemes.time import ExplicitEuler class AdvectionScheme(Upwind, ExplicitEuler): def CFL( self, cells: MeshCellSet, CFL_value: float, ) -> float: U_abs = np.linalg.norm(self.problem.V) dx = np.min(cells.surfaces) return CFL_value * dx / U_abs scheme = AdvectionScheme(AdvectionProblem()) ``` ```python from josie.solver import Solver solver = Solver(mesh, AdvectionState, scheme) ``` We need now to define an initialization function to initialize the domain ```python def init_fun(cells: MeshCellSet): xc = cells.centroids[..., 0] xc_r = np.where(xc >= 0.5) xc_l = np.where(xc < 0.5) cells.values[xc_r[0], xc_r[1], ...] = 1 cells.values[xc_l[0], xc_l[1], ...] = 0 ``` ### Run the simulation Let's initialize the solver state ```python solver.init(init_fun) # Plotting stuff import matplotlib.pyplot as plt import copy from matplotlib.animation import ArtistAnimation solution = [] t = 0 final_time = 1 CFL = 0.5 ``` Let's iterate in time choosing a time-based writing strategy (save every `$ 0.01s $`) to memory ```python from josie.io.write.strategy import TimeStrategy from josie.io.write.writer import MemoryWriter strategy = TimeStrategy(dt_save=0.01) writer = MemoryWriter(strategy, solver, final_time, CFL) writer.solve() ``` ```python fig, ax = plt.subplots() ims = [] for solution in writer.data: cells = solution.mesh.cells x = cells.centroids[..., 0] x = x.reshape(x.size) phi = cells.values[..., AdvectionFields.phi] (im,) = ax.plot(x, phi, 'ko-') ims.append([im]) ani = ArtistAnimation(fig, ims, interval=50) ``` ![png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAANbklEQVR4nO3df6jd9X3H8efLZK6MWR3LLZQkNZZFaHAD5SKOwurQjZg/kj+6lQSk6wiGdrMMWgYOhyvpX66sg0K2NmPiWqg27R/lQlMC6xRBGpcrWmsiltvUNjeVeWud/4jVsPf+OMdxdr0355vke8/J/eT5gMA53/PxnPcn5+bpyfmRk6pCkrT+XTXtASRJ/TDoktQIgy5JjTDoktQIgy5Jjdg4rRvetGlTbdu2bVo3L0nr0tNPP/2LqppZ6bKpBX3btm3Mz89P6+YlaV1K8tPVLvMpF0lqhEGXpEYYdElqhEGXpEYYdElqxNigJ3koyStJnl/l8iT5UpKFJM8luaX/MSVJ43R5hP4wsPM8l98FbB/+OgD886WPJUm6UGODXlVPAL88z5I9wFdr4DhwXZL39zWgJKmbPp5D3wycGTm/ODz2LkkOJJlPMr+0tNTDTUuS3jHRF0Wr6nBVzVbV7MzMip9clSRdpD6CfhbYOnJ+y/CYJGmC+gj6HPDx4btdbgNer6qXe7heSdIFGPuPcyV5BLgd2JRkEfg74NcAqurLwFFgF7AAvAH8+VoNK0la3digV9W+MZcX8Je9TSRJuih+UlSSGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGtEp6El2JnkxyUKS+1a4/ANJHkvyTJLnkuzqf1RJ0vmMDXqSDcAh4C5gB7AvyY5ly/4WOFJVNwN7gX/qe1BJ0vl1eYR+K7BQVaer6i3gUWDPsjUFvHd4+lrg5/2NKEnqokvQNwNnRs4vDo+N+hxwd5JF4Cjw6ZWuKMmBJPNJ5peWli5iXEnSavp6UXQf8HBVbQF2AV9L8q7rrqrDVTVbVbMzMzM93bQkCboF/SywdeT8luGxUfuBIwBV9X3gPcCmPgaUJHXTJegngO1JbkhyNYMXPeeWrfkZcAdAkg8xCLrPqUjSBI0NelWdA+4FjgEvMHg3y8kkB5PsHi77LHBPkh8AjwCfqKpaq6ElSe+2scuiqjrK4MXO0WMPjJw+BXy439EkSRfCT4pKUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1olPQk+xM8mKShST3rbLmY0lOJTmZ5Ov9jilJGmfjuAVJNgCHgD8CFoETSeaq6tTImu3A3wAfrqrXkrxvrQaWJK2syyP0W4GFqjpdVW8BjwJ7lq25BzhUVa8BVNUr/Y4pSRqnS9A3A2dGzi8Oj426EbgxyZNJjifZudIVJTmQZD7J/NLS0sVNLElaUV8vim4EtgO3A/uAf0ly3fJFVXW4qmaranZmZqanm5YkQbegnwW2jpzfMjw2ahGYq6q3q+onwI8YBF6SNCFdgn4C2J7khiRXA3uBuWVrvs3g0TlJNjF4CuZ0f2NKksYZG/SqOgfcCxwDXgCOVNXJJAeT7B4uOwa8muQU8Bjw11X16loNLUl6t1TVVG54dna25ufnp3LbkrReJXm6qmZXusxPikpSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIzoFPcnOJC8mWUhy33nWfTRJJZntb0RJUhdjg55kA3AIuAvYAexLsmOFddcAfwU81feQkqTxujxCvxVYqKrTVfUW8CiwZ4V1nwceBN7scT5JUkddgr4ZODNyfnF47P8kuQXYWlXfOd8VJTmQZD7J/NLS0gUPK0la3SW/KJrkKuCLwGfHra2qw1U1W1WzMzMzl3rTkqQRXYJ+Ftg6cn7L8Ng7rgFuAh5P8hJwGzDnC6OSNFldgn4C2J7khiRXA3uBuXcurKrXq2pTVW2rqm3AcWB3Vc2vycSSpBWNDXpVnQPuBY4BLwBHqupkkoNJdq/1gJKkbjZ2WVRVR4Gjy449sMra2y99LEnShfKTopLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY3oFPQkO5O8mGQhyX0rXP6ZJKeSPJfke0mu739USdL5jA16kg3AIeAuYAewL8mOZcueAWar6veAbwF/3/egkqTz6/II/VZgoapOV9VbwKPAntEFVfVYVb0xPHsc2NLvmJKkcboEfTNwZuT84vDYavYD313pgiQHkswnmV9aWuo+pSRprF5fFE1yNzALfGGly6vqcFXNVtXszMxMnzctSVe8jR3WnAW2jpzfMjz2/yS5E7gf+EhV/aqf8SRJXXV5hH4C2J7khiRXA3uBudEFSW4GvgLsrqpX+h9TkjTO2KBX1TngXuAY8AJwpKpOJjmYZPdw2ReA3wS+meTZJHOrXJ0kaY10ecqFqjoKHF127IGR03f2PJck6QL5SVFJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJakSnoCfZmeTFJAtJ7lvh8l9P8o3h5U8l2db7pJKk8xob9CQbgEPAXcAOYF+SHcuW7Qdeq6rfAf4ReLDvQSVJ59flEfqtwEJVna6qt4BHgT3L1uwB/m14+lvAHUnS35iSpHG6BH0zcGbk/OLw2Iprquoc8Drw28uvKMmBJPNJ5peWli5uYknSiib6omhVHa6q2aqanZmZmeRNS1LzugT9LLB15PyW4bEV1yTZCFwLvNrHgJKkbroE/QSwPckNSa4G9gJzy9bMAX82PP0nwH9UVfU3piRpnI3jFlTVuST3AseADcBDVXUyyUFgvqrmgH8FvpZkAfglg+hLkiZobNABquoocHTZsQdGTr8J/Gm/o0mSLoSfFJWkRhh0SWqEQZekRhh0SWpEpvXuwiRLwE8v8j/fBPyix3HWA/d8ZXDPV4ZL2fP1VbXiJzOnFvRLkWS+qmanPcckuecrg3u+MqzVnn3KRZIaYdAlqRHrNeiHpz3AFLjnK4N7vjKsyZ7X5XPokqR3W6+P0CVJyxh0SWrEZR30K/HLqTvs+TNJTiV5Lsn3klw/jTn7NG7PI+s+mqSSrPu3uHXZc5KPDe/rk0m+PukZ+9bhZ/sDSR5L8szw53vXNObsS5KHkryS5PlVLk+SLw1/P55Lcssl32hVXZa/GPxTvT8GPghcDfwA2LFszV8AXx6e3gt8Y9pzT2DPfwj8xvD0p66EPQ/XXQM8ARwHZqc99wTu5+3AM8BvDc+/b9pzT2DPh4FPDU/vAF6a9tyXuOc/AG4Bnl/l8l3Ad4EAtwFPXeptXs6P0K/EL6ceu+eqeqyq3hiePc7gG6TWsy73M8DngQeBNyc53Brpsud7gENV9RpAVb0y4Rn71mXPBbx3ePpa4OcTnK93VfUEg++HWM0e4Ks1cBy4Lsn7L+U2L+eg9/bl1OtIlz2P2s/g//Dr2dg9D/8qurWqvjPJwdZQl/v5RuDGJE8mOZ5k58SmWxtd9vw54O4kiwy+f+HTkxltai70z/tYnb7gQpefJHcDs8BHpj3LWkpyFfBF4BNTHmXSNjJ42uV2Bn8LeyLJ71bVf09zqDW2D3i4qv4hye8z+Ba0m6rqf6Y92HpxOT9CvxK/nLrLnklyJ3A/sLuqfjWh2dbKuD1fA9wEPJ7kJQbPNc6t8xdGu9zPi8BcVb1dVT8BfsQg8OtVlz3vB44AVNX3gfcw+EesWtXpz/uFuJyDfiV+OfXYPSe5GfgKg5iv9+dVYcyeq+r1qtpUVduqahuD1w12V9X8dMbtRZef7W8zeHROkk0MnoI5PcEZ+9Zlzz8D7gBI8iEGQV+a6JSTNQd8fPhul9uA16vq5Uu6xmm/EjzmVeJdDB6Z/Bi4f3jsIIM/0DC4w78JLAD/CXxw2jNPYM//DvwX8Ozw19y0Z17rPS9b+zjr/F0uHe/nMHiq6RTwQ2DvtGeewJ53AE8yeAfMs8AfT3vmS9zvI8DLwNsM/sa1H/gk8MmR+/jQ8Pfjh338XPvRf0lqxOX8lIsk6QIYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEb8L0OdxLwyHnilAAAAAElFTkSuQmCC )