{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Computação Científica com Julia\n",
"\n",
" * Computação científica\n",
" * Julia\n",
"## Como são as coisas no IPT hoje?\n",
"\n",
" * Planilha eletrônica\n",
" * Excel\n",
" * Deus nos acuda!\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Será que isso é o suficiente???\n",
"\n",
" * Como detectar erros em planilhas?\n",
" * Como reutilizar cálculos já feitos?\n",
" * Como fazer algo mais complexo?\n",
" * Como fazer **gráficos decentes**?\n",
"\n",
"## O problema do Excel é mais um problema de cultura do que qualquer outra coisa!"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# O que eu quero dizer com cultura:\n",
"\n",
" * Modelagem matemática de problemas se resume ao que se fazia 100 anos atrás.\n",
" * Pouco reaproveitamento \n",
" * Pouco compartilhamento de trabalho já realizado\n",
" * Endeusamento do que já funciona\n",
" * Pouca discussão e ossificação do conhecimento\n",
" \n",
"## Exemplo do que já ocorreu no passado\n",
"\n",
" * Túnel de vento (tuninho) circa 1990 - Mestrado do Nilson\n",
" * Algor - escoamento potencial\n",
" * STAN5 (ou 7 de acordo com o Nilson)\n",
" * Modelo em escala reduzida\n",
" * Medição com fio quente"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Minha opinião\n",
"\n",
" * Compreender o que ocorre numa bancada ou em um trabalho externo\n",
" * Dar ao computador o seu lugar de direito: um escravo muito burro mas muito forte!\n",
" * Muito do que nós fazemos é software. **Engenharia de Software**\n",
" * Conhecer o que se faz aqui dentro e fora\n",
" * Compartilhar este conhecimento"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Ideal\n",
"\n",
" * Todas as bancadas modeladas\n",
" * Modelos hidráulicos **e**\n",
" * Modelos CFD\n",
" * Conhecimento profundo do que ocorre de fato nas bancadas.\n",
" * Um modelo das plantas em campo **antes** de medir.\n",
" * Cada área ter um conjunto de ferramentas bem conhecidas \n",
" * Ser \"piloto\" de software não vai nos levar longe \n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Computação Científica\n",
"\n",
" * Métodos básicos \n",
" * Sistemas lineares\n",
" * Aproximação\n",
" * Equações diferenciais\n",
" * Sistemas não lineares\n",
" * Mais importante do que o básico, é saber como usar isso.\n",
"\n",
" * Métodos numérico rápidos e eficientes\n",
" * Muito trabalhoso\n",
" * Muito difícil\n",
" * Coisa de profissa. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Porque Julia\n",
"\n",
"1. Ambiente interativo\n",
" * Muito melhor que FORTRAN ou C/C++/Pascal/Java\n",
" * Parecido com R/Python/Matlab/Mathematica/Scilab\n",
"2. Desempenho\n",
" * Consegue desempenho de C/FORTRAN\n",
" * Mas não precisa vetorizar\n",
" * Praticamente tudo é implementado em Julia!\n",
"3. Biblioteca extensa\n",
" * Visualização, Álgebra Linear, equações diferenciais, etc\n",
" * Ainda não tem tudo o que Matlab tem\n",
" * Muito fácil chamar Python (que hoje tem um ecossistema impressionante)\n",
"4. A linguagem é muito boa!\n",
"\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Engenharia de Software\n",
"\n",
" * Controle de versões e repositórios de software (git e github)\n",
" * Documentação no código (*literate programming*)\n",
" * Testes unitários\n",
" * Pensar no ambiente de trabalho:\n",
" * Editor\n",
" * Interface Notebook - Jupyter\n",
" * REPL (terminal)\n",
" * *Muito do que nós fazemos é programação exploratória*"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Estrutura do curso\n",
"\n",
"Objetivo final:\n",
"$$\n",
"\\nabla\\cdot\\pmb{u} = 0\\\\\n",
"\\frac{\\partial \\pmb{u}}{\\partial t} + \\pmb{u}\\cdot\\nabla\\pmb{u} = \n",
"-\\nabla p + \\frac{1}{Re}\\nabla^2\\pmb{u}\n",
"$$\n",
"\n",
"Precisamos saber resolver este tipo de equação:\n",
"$$\n",
"\\nabla^2 \\phi = f\\\\\n",
"\\frac{\\partial \\phi}{\\partial t} = \\alpha\\nabla^2 \n",
"\\phi\\\\\n",
"\\pmb{u}\\cdot\\phi = \\alpha\\nabla^2\\phi\n",
"$$\n",
"\n",
"Vamos focar no problema 1D para depois tentar avançar em outras direções"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Em paralelo: use algum software para resolver um problema concreto:\n",
" \n",
" * [OpenFoam](https://www.openfoam.com/)\n",
" * [SU2](https://su2code.github.io/)\n",
" * [CalculiX](http://www.calculix.de/)\n",
" * [Code Saturne](https://www.code-saturne.org/)\n",
" * [Code Aster](https://www.code-aster.org/)\n",
" * [Fenics](https://fenicsproject.org/)\n",
" * [FreeFem++](https://freefem.org/)\n",
" * [Dedalus Project](http://dedalus-project.org/)\n",
" * [Nektar](https://www.nektar.info/)\n",
" * [MEEP](https://meep.readthedocs.io/en/latest/)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Roteiro para $\\nabla^2 u = f$\n",
"\n",
" * $u^\\delta(x) = \\sum \\hat{u}_k \\phi_k(x) \\approx u(x)$ - Aproximação e interpolação\n",
" * $\\nabla^2 u^\\delta - f = \\varepsilon(x)$ - Erro da equação diferencial\n",
" * Resíduos ponderados $\\int_\\Omega w_i(x) \\nabla^2 u^\\delta \\:dx = \\int_\\Omega f w_i(x)\\:dx$\n",
" * Formulação fraca: $ -\\int_\\Omega \\nabla w_i(x) \\cdot \\nabla u^\\delta \\:dx + \\int_{\\partial\\Omega} w_i\\frac{\\partial u}{\\partial n} = \\int_\\Omega f w_i(x)\\:dx$\n",
" * Galerkin: $w_i(x) = \\phi_i(x)$\n",
" * Sistema linear: $\\left[A\\right]\\cdot\\left\\{\\hat{u}\\right\\} = \\left\\{f\\right\\}$, \n",
" $$A_{i,k} = -\\int_\\Omega\\nabla\\phi_i\\cdot\\nabla\\phi_k\\:dx + \\int_{\\partial\\Omega}\\phi_i\\frac{\\partial u}{\\partial n}\\\\\n",
" f_i = \\int_\\Omega f\\phi_i\\:dx\n",
" $$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Problema 1D\n",
"\n",
" 1. Interpolação e aproximação: polinomial, senos e cossenos\n",
" 2. Derivadas: simbólica, numérica, diferenciação automática\n",
" 3. Quadratura, Simbólica, Numérica\n",
" 4. Sistemas lineares: Métodos diretos, métodos iterativos, *Explorar a estrutura da matriz*\n",
" 5. Equações diferenciais ordinárias: problemas de valor inicial\n",
" 6. Sistemas de equações não lineares: Newton-Raphson e outros\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Introduzindo Julia\n",
"\n",
" * http://julialang.org\n",
" * Lista de discussão https://discourse.julialang.org/\n",
" * Canal do youtube https://www.youtube.com/user/JuliaLanguage\n",
" * Documentação: https://docs.julialang.org"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Caluladora\n",
"1+1\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"21"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function soma(a,b) \n",
" return a+b\n",
"end\n",
"soma2(a,b) = a+2b\n",
"\n",
"x = 1\n",
"y = 10\n",
"soma2(x,y)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Porque Julia consegue ser rápido como C?"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\t.text\n",
"; ┌ @ In[10]:4 within `soma2'\n",
"; │┌ @ In[10]:4 within `+'\n",
"\tleaq\t(%rdi,%rsi,2), %rax\n",
"; │└\n",
"\tretq\n",
"\tnopw\t%cs:(%rax,%rax)\n",
"\tnop\n",
"; └\n"
]
}
],
"source": [
"code_native(soma2, (Int,Int))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\t.text\n",
"; ┌ @ In[10]:4 within `soma2'\n",
"; │┌ @ promotion.jl:312 within `*' @ In[10]:4\n",
"\tvaddsd\t%xmm1, %xmm1, %xmm1\n",
"; │└\n",
"; │┌ @ float.jl:401 within `+'\n",
"\tvaddsd\t%xmm0, %xmm1, %xmm0\n",
"; │└\n",
"\tretq\n",
"\tnopl\t(%rax)\n",
"; └\n"
]
}
],
"source": [
"code_native(soma2, (Float64,Float64))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"df1b (generic function with 1 method)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using ForwardDiff\n",
"\n",
"function f1(x)\n",
" \n",
" res = one(x)\n",
" for i = 1:5\n",
" res = x + res*x\n",
" end\n",
" return res\n",
"end\n",
"df1 = x -> ForwardDiff.derivative(f1, x)\n",
"f1b(x) = x + x*(one(x) + x*(one(x) + x * (one(x)+ 2x)))\n",
"df1b(x) = one(x) + 2x + 3x^2 + 4x^3 + 10x^4\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df1(0.5) - df1b(0.5)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.736 ns (0 allocations: 0 bytes)\n",
" 0.024 ns (0 allocations: 0 bytes)\n",
" 14.957 ns (0 allocations: 0 bytes)\n",
" 0.024 ns (0 allocations: 0 bytes)\n"
]
},
{
"data": {
"text/plain": [
"20"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using BenchmarkTools\n",
"\n",
"x - 0.5\n",
"@btime f1($x)\n",
"@btime f1b($x)\n",
"@btime df1($x)\n",
"@btime df1b($x)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Sistema de equações diferenciais ordinárias:\n",
"\n",
"Equações de Lorenz\n",
"$$\n",
"\\begin{align}\n",
"\\frac{dx}{dt} &= σ(y-x) \\\\\n",
"\\frac{dy}{dt} &= x(ρ-z) - y \\\\\n",
"\\frac{dz}{dt} &= xy - βz \\\\\n",
"\\end{align}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Plots.GRBackend()"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using DifferentialEquations\n",
"using Plots\n",
"gr()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"lorenz (generic function with 1 method)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function lorenz(du,u,p,t)\n",
" du[1] = 10.0*(u[2]-u[1])\n",
" du[2] = u[1]*(28.0-u[3]) - u[2]\n",
" du[3] = u[1]*u[2] - (8/3)*u[3]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"u0 = [1.0;0.0;0.0]\n",
"tspan = (0.0,100.0)\n",
"prob = ODEProblem(lorenz,u0,tspan)\n",
"sol = solve(prob);"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gr()\n",
"plot(sol,vars=(1,2,3))"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"using Psychro\n",
"using Unitful"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"383282.6011056348 cm^3 lb^-1\n",
"17.134466939662275 °C\n",
"0.755640835595789\n"
]
}
],
"source": [
"println(volume(MoistAir, 20.0u\"°C\", DewPoint, 60.0u\"°F\", 1.0u\"atm\", u\"cm^3/lb\"))\n",
"println(wetbulb(MoistAir, 20.0u\"°C\", DewPoint, 60.0u\"°F\", 1.0u\"atm\"))\n",
"println(relhum(MoistAir, 20.0u\"°C\", DewPoint, 60.0u\"°F\", 1.0u\"atm\"))\n"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using ACME\n",
"\n",
"circ = @circuit begin\n",
" j_in = voltagesource(), [-] ⟷ gnd\n",
" r1 = resistor(1e3), [1] ⟷ j_in[+]\n",
" c1 = capacitor(47e-9), [1] ⟷ r1[2], [2] ⟷ gnd\n",
" d1 = diode(is=1e-15), [+] ⟷ r1[2], [-] ⟷ gnd\n",
" d2 = diode(is=1.8e-15), [+] ⟷ gnd, [-] ⟷ r1[2]\n",
" j_out = voltageprobe(), [+] ⟷ r1[2], [-] ⟷ gnd\n",
"end\n",
"model = DiscreteModel(circ, 1/4410)\n",
"t = 1/44100 * (0:44099)'\n",
"y = run!(model, sin.(2π*100 .* t));\n",
"plot(t[1,:],y[1,:])"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"using PyCall"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.8414709848078965\n"
]
},
{
"data": {
"text/plain": [
"3×4 Array{Float64,2}:\n",
" 0.402438 0.00215392 1.24917 -0.291974\n",
" 1.63742 -0.49273 -0.0018262 0.522377\n",
" 1.17184 0.505693 0.11604 -1.42486"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"math = pyimport(\"math\")\n",
"println(math.sin(1))\n",
"\n",
"np = pyimport(\"numpy\")\n",
"nprand = pyimport(\"numpy.random\")\n",
"nprand.randn(3,4)\n"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Fun(Chebyshev(0..10),[-0.05037748828223824, -0.6282471738193445, -0.02641752529939319, -0.6996910463813416, 0.22195565278327428, 0.5008280175468399, -0.07434713931875803, -0.10236787078571936, 0.01044172796045064, 0.010587067010913753 … -4.1863293741963567e-10, 1.5716757540826944e-11, 6.412940310800903e-12, -2.1829925116784714e-13, -8.114196016050454e-14, 2.5268751058932945e-15, 8.625777787526683e-16, -2.4745528200200695e-17, -7.548648766675401e-18, -2.0553989296917155e-19])"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra, SpecialFunctions, Plots, ApproxFun\n",
"x = Fun(identity,0..10)\n",
"f = sin(x^2)\n",
"g = cos(x)"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h = f + g^2\n",
"r = roots(h)\n",
"rp = roots(h')\n",
"\n",
"using Plots\n",
"plot(h)\n",
"scatter!(r,h.(r))\n",
"scatter!(rp,h.(rp))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 1.4.1",
"language": "julia",
"name": "julia-1.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.4.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}