Commit d7b3e0e8 authored by Jean-Noël Grad's avatar Jean-Noël Grad
Browse files

Fix typos in comments and docstrings

parent 85cb8eeb
...@@ -1452,7 +1452,7 @@ int main( int argc, char **argv ) ...@@ -1452,7 +1452,7 @@ int main( int argc, char **argv )
WALBERLA_LOG_INFO_ON_ROOT("Refreshing blockforest...") WALBERLA_LOG_INFO_ON_ROOT("Refreshing blockforest...")
// check refinement criterions and refine/coarsen if necessary // check refinement criteria and refine/coarsen if necessary
uint_t stampBefore = blocks->getBlockForest().getModificationStamp(); uint_t stampBefore = blocks->getBlockForest().getModificationStamp();
blocks->refresh(); blocks->refresh();
uint_t stampAfter = blocks->getBlockForest().getModificationStamp(); uint_t stampAfter = blocks->getBlockForest().getModificationStamp();
...@@ -2090,7 +2090,7 @@ int main( int argc, char **argv ) ...@@ -2090,7 +2090,7 @@ int main( int argc, char **argv )
WALBERLA_LOG_INFO_ON_ROOT("Refreshing blockforest...") WALBERLA_LOG_INFO_ON_ROOT("Refreshing blockforest...")
// check refinement criterions and refine/coarsen if necessary // check refinement criteria and refine/coarsen if necessary
uint_t stampBefore = blocks->getBlockForest().getModificationStamp(); uint_t stampBefore = blocks->getBlockForest().getModificationStamp();
blocks->refresh(); blocks->refresh();
uint_t stampAfter = blocks->getBlockForest().getModificationStamp(); uint_t stampAfter = blocks->getBlockForest().getModificationStamp();
......
...@@ -929,7 +929,7 @@ int main( int argc, char **argv ) ...@@ -929,7 +929,7 @@ int main( int argc, char **argv )
if( !useStaticRefinement && refinementCheckFrequency == 0 && numberOfLevels != 1 ) if( !useStaticRefinement && refinementCheckFrequency == 0 && numberOfLevels != 1 )
{ {
// determine check frequency automatically based on maximum admissable velocity and block sizes // determine check frequency automatically based on maximum admissible velocity and block sizes
real_t uMax = real_t(0.1); real_t uMax = real_t(0.1);
real_t refinementCheckFrequencyFinestLevel = ( overlap + real_c(blockSize) - real_t(2) * real_t(FieldGhostLayers) * dx) / uMax; real_t refinementCheckFrequencyFinestLevel = ( overlap + real_c(blockSize) - real_t(2) * real_t(FieldGhostLayers) * dx) / uMax;
refinementCheckFrequency = uint_c( refinementCheckFrequencyFinestLevel / real_t(lbmTimeStepsPerTimeLoopIteration)); refinementCheckFrequency = uint_c( refinementCheckFrequencyFinestLevel / real_t(lbmTimeStepsPerTimeLoopIteration));
...@@ -1252,7 +1252,7 @@ int main( int argc, char **argv ) ...@@ -1252,7 +1252,7 @@ int main( int argc, char **argv )
(*velocityCommunicationScheme)(); (*velocityCommunicationScheme)();
} }
// check refinement criterions and refine/coarsen if necessary // check refinement criteria and refine/coarsen if necessary
uint_t stampBefore = blocks->getBlockForest().getModificationStamp(); uint_t stampBefore = blocks->getBlockForest().getModificationStamp();
blocks->refresh(); blocks->refresh();
uint_t stampAfter = blocks->getBlockForest().getModificationStamp(); uint_t stampAfter = blocks->getBlockForest().getModificationStamp();
......
...@@ -6,7 +6,7 @@ import os ...@@ -6,7 +6,7 @@ import os
class Parameter: class Parameter:
def __init__(self, name, type, defValue="", comment=""): def __init__(self, name, type, defValue="", comment=""):
"""Propery of a data strcuture """Property of a data structure
Parameters Parameters
---------- ----------
......
...@@ -878,7 +878,7 @@ int main( int argc, char **argv ) ...@@ -878,7 +878,7 @@ int main( int argc, char **argv )
real_t defaultOmegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, real_t(1)); real_t defaultOmegaBulk = lbm_mesapd_coupling::omegaBulkFromOmega(omega, real_t(1));
shared_ptr<OmegaBulkAdapter_T> omegaBulkAdapter = make_shared<OmegaBulkAdapter_T>(blocks, omegaBulkFieldID, accessor, defaultOmegaBulk, omegaBulk, adaptionLayerSize, sphereSelector); shared_ptr<OmegaBulkAdapter_T> omegaBulkAdapter = make_shared<OmegaBulkAdapter_T>(blocks, omegaBulkFieldID, accessor, defaultOmegaBulk, omegaBulk, adaptionLayerSize, sphereSelector);
timeloopAfterParticles.add() << Sweep( makeSharedSweep(omegaBulkAdapter), "Omega Bulk Adapter"); timeloopAfterParticles.add() << Sweep( makeSharedSweep(omegaBulkAdapter), "Omega Bulk Adapter");
// initally adapt // initially adapt
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) { for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
(*omegaBulkAdapter)(blockIt.get()); (*omegaBulkAdapter)(blockIt.get());
} }
......
...@@ -843,7 +843,7 @@ int main( int argc, char **argv ) ...@@ -843,7 +843,7 @@ int main( int argc, char **argv )
auto sphereShape = ss->create<mesa_pd::data::Sphere>( diameter * real_t(0.5) ); auto sphereShape = ss->create<mesa_pd::data::Sphere>( diameter * real_t(0.5) );
ss->shapes[sphereShape]->updateMassAndInertia(densityRatio); ss->shapes[sphereShape]->updateMassAndInertia(densityRatio);
std::mt19937 randomGenerator (static_cast<unsigned int>(2610)); // fixed seed: quasi-random and reproducable std::mt19937 randomGenerator (static_cast<unsigned int>(2610)); // fixed seed: quasi-random and reproducible
for( uint_t nSed = 0; nSed < numberOfSediments; ++nSed ) for( uint_t nSed = 0; nSed < numberOfSediments; ++nSed )
{ {
...@@ -962,7 +962,7 @@ int main( int argc, char **argv ) ...@@ -962,7 +962,7 @@ int main( int argc, char **argv )
if(currentPhase == 1) if(currentPhase == 1)
{ {
// damp velocites to avoid too large ones // damp velocities to avoid too large ones
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
[](const size_t idx, ParticleAccessor_T& ac){ [](const size_t idx, ParticleAccessor_T& ac){
ac.setLinearVelocity(idx, ac.getLinearVelocity(idx) * real_t(0.5)); ac.setLinearVelocity(idx, ac.getLinearVelocity(idx) * real_t(0.5));
......
...@@ -573,7 +573,7 @@ int main( int argc, char **argv ) ...@@ -573,7 +573,7 @@ int main( int argc, char **argv )
if(maxPenetrationDepth < overlapLimit) break; if(maxPenetrationDepth < overlapLimit) break;
// reset velocites to avoid too large ones // reset velocities to avoid too large ones
ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor, ps->forEachParticle( useOpenMP, mesa_pd::kernel::SelectLocal(), *accessor,
[](const size_t idx, ParticleAccessor_T& ac){ [](const size_t idx, ParticleAccessor_T& ac){
......
...@@ -6,7 +6,7 @@ import os ...@@ -6,7 +6,7 @@ import os
class Parameter: class Parameter:
def __init__(self, name, type, defValue=""): def __init__(self, name, type, defValue=""):
"""Propery of a data strcuture """Property of a data structure
Parameters Parameters
---------- ----------
......
...@@ -1064,7 +1064,7 @@ void keepInflowOutflowAtTheSameLevel( std::vector< std::pair< const Block *, uin ...@@ -1064,7 +1064,7 @@ void keepInflowOutflowAtTheSameLevel( std::vector< std::pair< const Block *, uin
uint_t maxInflowLevel( uint_t(0) ); uint_t maxInflowLevel( uint_t(0) );
uint_t maxOutflowLevel( uint_t(0) ); uint_t maxOutflowLevel( uint_t(0) );
// In addtion to keeping in- and outflow blocks at the same level, this callback also // In addition to keeping in- and outflow blocks at the same level, this callback also
// prevents these blocks from coarsening. // prevents these blocks from coarsening.
for( auto it = minTargetLevels.begin(); it != minTargetLevels.end(); ++it ) for( auto it = minTargetLevels.begin(); it != minTargetLevels.end(); ++it )
...@@ -2569,14 +2569,14 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod ...@@ -2569,14 +2569,14 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
blockforest::DynamicDiffusionBalance< blockforest::NoPhantomData >( maxIterations, flowIterations ) ); blockforest::DynamicDiffusionBalance< blockforest::NoPhantomData >( maxIterations, flowIterations ) );
} }
// add callback functions which are executed after all block data was unpakced after the dynamic load balancing // add callback functions which are executed after all block data was unpacked after the dynamic load balancing
// for blocks that have *not* migrated: store current flag field state (required for lbm::PostProcessing) // for blocks that have *not* migrated: store current flag field state (required for lbm::PostProcessing)
blockforest.addRefreshCallbackFunctionAfterBlockDataIsUnpacked( lbm::MarkerFieldGenerator< LatticeModel_T, field::FlagFieldEvaluationFilter<FlagField_T> >( blockforest.addRefreshCallbackFunctionAfterBlockDataIsUnpacked( lbm::MarkerFieldGenerator< LatticeModel_T, field::FlagFieldEvaluationFilter<FlagField_T> >(
pdfFieldId, markerDataId, flagFieldFilter ) ); pdfFieldId, markerDataId, flagFieldFilter ) );
// (re)set boundaries = (re)initialize flag field for every block with respect to the new block structure (the size of neighbor blocks might have changed) // (re)set boundaries = (re)initialize flag field for every block with respect to the new block structure (the size of neighbor blocks might have changed)
blockforest.addRefreshCallbackFunctionAfterBlockDataIsUnpacked( blockforest::BlockForest::RefreshCallbackWrappper( boundarySetter ) ); blockforest.addRefreshCallbackFunctionAfterBlockDataIsUnpacked( blockforest::BlockForest::RefreshCallbackWrappper( boundarySetter ) );
// treat boundary-fluid cell convertions // treat boundary-fluid cell conversions
blockforest.addRefreshCallbackFunctionAfterBlockDataIsUnpacked( lbm::PostProcessing< LatticeModel_T, field::FlagFieldEvaluationFilter<FlagField_T> >( blockforest.addRefreshCallbackFunctionAfterBlockDataIsUnpacked( lbm::PostProcessing< LatticeModel_T, field::FlagFieldEvaluationFilter<FlagField_T> >(
pdfFieldId, markerDataId, flagFieldFilter ) ); pdfFieldId, markerDataId, flagFieldFilter ) );
// (re)set velocity field (velocity field data is not migrated!) // (re)set velocity field (velocity field data is not migrated!)
......
...@@ -295,7 +295,7 @@ class DummySweep ...@@ -295,7 +295,7 @@ class DummySweep
void emptyFunction() {} void emptyFunction() {}
//******************************************************************************************************************* //*******************************************************************************************************************
/*!\brief Simualtion of a strongly heterogeneous sized particulate flow system using combined resolved and unresolved /*!\brief Simulation of a strongly heterogeneous sized particulate flow system using combined resolved and unresolved
* methods. * methods.
* *
* For the coupling of resolved particles the Momentum Exchange Method (MEM) is used, whereas for the * For the coupling of resolved particles the Momentum Exchange Method (MEM) is used, whereas for the
......
...@@ -599,7 +599,7 @@ int main(int argc, char** argv) { ...@@ -599,7 +599,7 @@ int main(int argc, char** argv) {
WALBERLA_CHECK(!(useCurlCriterion && useVorticityCriterion), WALBERLA_CHECK(!(useCurlCriterion && useVorticityCriterion),
"Using curl and vorticity criterion together makes no sense."); "Using curl and vorticity criterion together makes no sense.");
// create base dir if it doesnt already exist // create base dir if it doesn't already exist
filesystem::path bpath(baseFolder); filesystem::path bpath(baseFolder);
if (!filesystem::exists(bpath)) { if (!filesystem::exists(bpath)) {
filesystem::create_directory(bpath); filesystem::create_directory(bpath);
......
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from pystencils.session import * from pystencils.session import *
import sympy as sp import sympy as sp
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
sp.init_printing() sp.init_printing()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Code Generation for the Heat Equation ### Code Generation for the Heat Equation
The <a target="_blank" href="https://en.wikipedia.org/wiki/Heat_equation">heat equation</a> which is a simple partial differential equation describing the flow of heat through a homogenous medium. We can write it as The <a target="_blank" href="https://en.wikipedia.org/wiki/Heat_equation">heat equation</a> which is a simple partial differential equation describing the flow of heat through a homogenous medium. We can write it as
$$ $$
\frac{\partial u}{\partial t} = \frac{\partial u}{\partial t} =
\kappa \left( \kappa \left(
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial x^2} +
\frac{\partial^2 u}{\partial y^2} \frac{\partial^2 u}{\partial y^2}
\right) \right)
$$ $$
where $\kappa$ is the medium's diffusion coefficient and $u(x, y, t)$ is the unknown temperature distribution at the coordinate $(x,y)$ at time $t$. where $\kappa$ is the medium's diffusion coefficient and $u(x, y, t)$ is the unknown temperature distribution at the coordinate $(x,y)$ at time $t$.
To discretize this equation using pystencils, we first need to define all the fields and other symbols involved. To discretize this equation using pystencils, we first need to define all the fields and other symbols involved.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
u, u_tmp = ps.fields("u, u_tmp: [2D]", layout='fzyx') u, u_tmp = ps.fields("u, u_tmp: [2D]", layout='fzyx')
kappa = sp.Symbol("kappa") kappa = sp.Symbol("kappa")
dx = sp.Symbol("dx") dx = sp.Symbol("dx")
dt = sp.Symbol("dt") dt = sp.Symbol("dt")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
We define the PDE using the pystencils building blocks for transient and spatial derivatives. The definition is implicitly equalled to zero. We use `ps.fd.transient` for a first derivative by time and `ps.fd.diff` to express the second derivatives. `ps.fd.diff` takes a field and a list of spatial dimensions in which the field should be differentiated. We define the PDE using the pystencils building blocks for transient and spatial derivatives. The definition is implicitly equalled to zero. We use `ps.fd.transient` for a first derivative by time and `ps.fd.diff` to express the second derivatives. `ps.fd.diff` takes a field and a list of spatial dimensions in which the field should be differentiated.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
heat_pde = ps.fd.transient(u) - kappa * ( ps.fd.diff( u, 0, 0 ) + ps.fd.diff( u, 1, 1 ) ) heat_pde = ps.fd.transient(u) - kappa * ( ps.fd.diff( u, 0, 0 ) + ps.fd.diff( u, 1, 1 ) )
heat_pde heat_pde
``` ```
%%%% Output: execute_result %%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAUAAAAAZCAYAAABAZvchAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJNUlEQVR4Ae2c/3XUOBDHl7wUkLur4EIH/KjgQgdwVwGkA3j8Bf/lhQ6ACiB0EKiAHx1wHcDRQe770Xr8ZK/slbxSVuvsvKeVJY1GoxnNaCzbe+vq6mpRI7x8+fKJ+PpX+cca+dvztJdATRKQnbxWOq2Jp13g5aBGJqXIh+Lr7t751aidPU+VSuBc9nJZKW/VsnWrtghQSjyStD4pvzsktQbnRO3HSh9V/jaEG1Ofm17MmDcdp4TMS9DcJT1p/tw1HSl/VRPfNeulxgjwQsp7PaRACRPH91wJp/dG6Z7q6DMJctObxMQN61RC5iVo7ppaJAPs4VQ5QUQVUL1eiABrSS9evDhW+m+IH7UdKZ3321X3VOmiX7+unJveuvH27VeLEjIvQXNXdSVZTLKFEvPdBb3UFgE+07bFLjYERH5PAjscfR6q/s5Qx4H63PQGhtlXexIoIfMSND2Wd+rSbKGGKLB6vdTmADnDGDvI/a7293J0v/wl2ZSp4/Y4BXLTSxl7J3El6xMl9DQVSsi8BM2p89tqv8YWOB76e1NGKtX1ptPq9D/slLZYkLAtevsyxIZwxqLDn+p3f6hvqD43vdAYM6wjspgcXZSQeQmaO643bOiR0pi9xEyxOl3HMJ2Ck80BahHy9IkobCoQvfHe3ygNtRN98ITY8D6r7oPKx/RXSoLc9JIGv6HIJWReguYOq4eIeOMIMMf8a9dLNgcoYfHqyjsl9wheObuHOwNQfmb1I0IlejOnFkQTDW6PcZLtC5+69m/JkhxgTnqixXxxzreVvqrc2X1V5kn1Y+WjcxTOrEHzz6pDhJWLpujMRYfcAhOQbBqUbLQWc+iltE6ynAE2THILi+AXKnPNDvRaiXA85msOFh+3sUEQTXc2qLx1fiCqDG2r+xrsHKjMTU9DPBdNnD98nvtDqh7HyEOaoPNTPW1EsEVAtJ8WIZxIVHxk1SHDZ6Y5Fx2aHRVbU+tUn1Evk3Wyjkfac0WALtzWpHkpma84fim3COhBDCPC+V0pGMGJFg6EW+Shl6NxnsB7foRvDuiHin+ozNPlFgrQw+F/bgZgvrYAbUzq3OZgFZY3vBDVtnPX9Sj/1tfPx/qo7ZXSVj+V0vipOkSnrCveawvqPSdN0ZqTDm2jxaauHVL1AoPqQ/DySLlvB5N1YpMWPTZ/7i6xSZPLma65Oz07FAIL7ZMSeSzAqG/QGDhOjx3eor5YWobH+H3HYW04MOj7Y1obOTvdB7WDg6H9UG634ie6vlTyHXFuejgw4w2jRcA+4Lz7dQv1Yc4PlHNg7UDXMfwbekqfC9Hmc6nOZtAh5BWEhx7huw/OqNR+2m9Q+Zvq27n02lNkzsK3scfWZU6ac9LhkB31VLIsblPXcKDx0fEd5a3zW3LmjruS7aqhyfphDT8TXecLGpqMh5/iA4pnOMBfKgR3WOsQkTMY0RfnXDgjzgOZTMqZF3ys7Fiig3Ccg1O+Amo3Q2GyANHTX+5KP2onKsUBHiuxyLPSYxzRhHdyeIG+Rb/UwTt1oWMAeDW+delglH9D6uVr+4gP5IADPFJy/PZodIrCCTm4heqJ8JFlZ1F1OvcKwmX+0ToUPoseZ8pYQchNU/ScTJTPQYdmR1GOUHPemq4b5SJzc3StvqfqRP1YN2+V/jQaLdHlBTbnnO1BryG5qAHMwIkwMLI3SjhUJkU040B1GMEUMGXaLWafBlEAxsLYxkt/J2FxE1UAuektqS5/iX5wss6YmgZnxKpbUbDaT+C7wVvoOoZ/Q3d5Yh/GQi/XDdEyT2CsBE2Gn4MOzdb8dZgg2o1Qo/XC2lVi8yYtuFa6Exg9Wifqjw2te+CIf3Dn0QeBwVKrnEFp4NaQ+wTUhiM0wfSbrQxTpjirW6gv9UFFqg3nwvh227XSvyHETujGz02voW8Zwu87Om69nWw09hMlx6NycPvzWsu/DeTlKX3YRPyjAI9MuctEmUcxUoJmM/AcdGi2xrq/VkjRC7hK7hhDTPKmCLerffuB/2idCBdnCl1ejQsCY1j7QRAjrbI18EA3i8Rua0C7DqC5Kozz3kAjQuqE6aLHTkGYy99mraPNgvAdRW56Iu+gw0fDI3MypSIHc3ooNXaB9vlfjjb+G+rD2Iy7DdhU5iGeS9Ccgw7RMefhttZCsitZl6oXbHkwgFJbik4IiMZodeZ92ClNK+BY+udYUHqs9I+UcH+gHRwfcBKcT62cUamO2+qfSoS2vtPo3+MPKRweWyHmpudNAsW/FX12IZ5AMybHAYT2PI16p2QAT31++2Uft+XfKps8pY8vux6ZssUEmUczUoKmBp+DDllz0U4gWuCRiCl6ES4Oa52zTtEJXPIi+CBoTF47cxHixg5QhIK3VM0Ag2Fonzvhc4aHMRMxrSgvhp5w7PyNHRCH6kOnnJseA4km/NvtuD92qA6HduQjpfBv/RL7IJchR2oki+XilfUQvSZiGMlNU/TmoENsKBSUxIg0C06CXtroT32wB84F+7aaohN8B0FXEESbQKRdgwdBrO1V8iQ56FATWDoTLruKA02Ya16RmWr4uekZXygZh9SH0fE0DyJkd2jsdRzt4+Ex3uju6OEOXbIYSTcepIfqdMj6kGJwKtjSpnAduubLKY6/gOdKU+3UEdAPd548XGx9gDWoDufH+V87RlX/CC3GUByv0PxmTE/J1d8cBLehCJjD1clGm5uezUl0eRLFi76tQmgbG09tTkZC4xWjdicb6wNNQDjB8Zat9fyKTxz1qRKLmPnyyk37bqeuk6EETZgYkqnqB9eg2orpULR54Nh5tzRZWNfYoZEFuuZF6C8qd6K/KayIBpuAyR+7xw8AHKV1/UCJP0LchKb+RPFC6eEmNHalr+Z5orTyB68x/KfKSPj8mexlDO09TvyfBEumVelQ/HxXOt7rME6HB0vHWNUvISyh8OxBuxHnFdzSEvFEg/CJjFbOSdcQYEdkp91DRgnUpEPxQvTHJ4+dO4qM050dqeocoJRHiMo7QRbCzk7o/oQ0T5wSt8GE7bHAQXE3lB/pKVwMgxfV94YxIqepTZLr1nUoHthE+UQ1+gudqfOdU7/qHCDClRI52/qufOUgc07Ct7lonjzmjwbht5/aRXbiX7RTI8ZI0ns0JFCBDllDobcN9goakcD/XTRMxcq6e8AAAAAASUVORK5CYII=) ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAUAAAAAZCAYAAABAZvchAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJNUlEQVR4Ae2c/3XUOBDHl7wUkLur4EIH/KjgQgdwVwGkA3j8Bf/lhQ6ACiB0EKiAHx1wHcDRQe770Xr8ZK/slbxSVuvsvKeVJY1GoxnNaCzbe+vq6mpRI7x8+fKJ+PpX+cca+dvztJdATRKQnbxWOq2Jp13g5aBGJqXIh+Lr7t751aidPU+VSuBc9nJZKW/VsnWrtghQSjyStD4pvzsktQbnRO3HSh9V/jaEG1Ofm17MmDcdp4TMS9DcJT1p/tw1HSl/VRPfNeulxgjwQsp7PaRACRPH91wJp/dG6Z7q6DMJctObxMQN61RC5iVo7ppaJAPs4VQ5QUQVUL1eiABrSS9evDhW+m+IH7UdKZ3321X3VOmiX7+unJveuvH27VeLEjIvQXNXdSVZTLKFEvPdBb3UFgE+07bFLjYERH5PAjscfR6q/s5Qx4H63PQGhtlXexIoIfMSND2Wd+rSbKGGKLB6vdTmADnDGDvI/a7293J0v/wl2ZSp4/Y4BXLTSxl7J3El6xMl9DQVSsi8BM2p89tqv8YWOB76e1NGKtX1ptPq9D/slLZYkLAtevsyxIZwxqLDn+p3f6hvqD43vdAYM6wjspgcXZSQeQmaO643bOiR0pi9xEyxOl3HMJ2Ck80BahHy9IkobCoQvfHe3ygNtRN98ITY8D6r7oPKx/RXSoLc9JIGv6HIJWReguYOq4eIeOMIMMf8a9dLNgcoYfHqyjsl9wheObuHOwNQfmb1I0IlejOnFkQTDW6PcZLtC5+69m/JkhxgTnqixXxxzreVvqrc2X1V5kn1Y+WjcxTOrEHzz6pDhJWLpujMRYfcAhOQbBqUbLQWc+iltE6ynAE2THILi+AXKnPNDvRaiXA85msOFh+3sUEQTXc2qLx1fiCqDG2r+xrsHKjMTU9DPBdNnD98nvtDqh7HyEOaoPNTPW1EsEVAtJ8WIZxIVHxk1SHDZ6Y5Fx2aHRVbU+tUn1Evk3Wyjkfac0WALtzWpHkpma84fim3COhBDCPC+V0pGMGJFg6EW+Shl6NxnsB7foRvDuiHin+ozNPlFgrQw+F/bgZgvrYAbUzq3OZgFZY3vBDVtnPX9Sj/1tfPx/qo7ZXSVj+V0vipOkSnrCveawvqPSdN0ZqTDm2jxaauHVL1AoPqQ/DySLlvB5N1YpMWPTZ/7i6xSZPLma65Oz07FAIL7ZMSeSzAqG/QGDhOjx3eor5YWobH+H3HYW04MOj7Y1obOTvdB7WDg6H9UG634ie6vlTyHXFuejgw4w2jRcA+4Lz7dQv1Yc4PlHNg7UDXMfwbekqfC9Hmc6nOZtAh5BWEhx7huw/OqNR+2m9Q+Zvq27n02lNkzsK3scfWZU6ac9LhkB31VLIsblPXcKDx0fEd5a3zW3LmjruS7aqhyfphDT8TXecLGpqMh5/iA4pnOMBfKgR3WOsQkTMY0RfnXDgjzgOZTMqZF3ys7Fiig3Ccg1O+Amo3Q2GyANHTX+5KP2onKsUBHiuxyLPSYxzRhHdyeIG+Rb/UwTt1oWMAeDW+delglH9D6uVr+4gP5IADPFJy/PZodIrCCTm4heqJ8JFlZ1F1OvcKwmX+0ToUPoseZ8pYQchNU/ScTJTPQYdmR1GOUHPemq4b5SJzc3StvqfqRP1YN2+V/jQaLdHlBTbnnO1BryG5qAHMwIkwMLI3SjhUJkU040B1GMEUMGXaLWafBlEAxsLYxkt/J2FxE1UAuektqS5/iX5wss6YmgZnxKpbUbDaT+C7wVvoOoZ/Q3d5Yh/GQi/XDdEyT2CsBE2Gn4MOzdb8dZgg2o1Qo/XC2lVi8yYtuFa6Exg9Wifqjw2te+CIf3Dn0QeBwVKrnEFp4NaQ+wTUhiM0wfSbrQxTpjirW6gv9UFFqg3nwvh227XSvyHETujGz02voW8Zwu87Om69nWw09hMlx6NycPvzWsu/DeTlKX3YRPyjAI9MuctEmUcxUoJmM/AcdGi2xrq/VkjRC7hK7hhDTPKmCLerffuB/2idCBdnCl1ejQsCY1j7QRAjrbI18EA3i8Rua0C7DqC5Kozz3kAjQuqE6aLHTkGYy99mraPNgvAdRW56Iu+gw0fDI3MypSIHc3ooNXaB9vlfjjb+G+rD2Iy7DdhU5iGeS9Ccgw7RMefhttZCsitZl6oXbHkwgFJbik4IiMZodeZ92ClNK+BY+udYUHqs9I+UcH+gHRwfcBKcT62cUamO2+qfSoS2vtPo3+MPKRweWyHmpudNAsW/FX12IZ5AMybHAYT2PI16p2QAT31++2Uft+XfKps8pY8vux6ZssUEmUczUoKmBp+DDllz0U4gWuCRiCl6ES4Oa52zTtEJXPIi+CBoTF47cxHixg5QhIK3VM0Ag2Fonzvhc4aHMRMxrSgvhp5w7PyNHRCH6kOnnJseA4km/NvtuD92qA6HduQjpfBv/RL7IJchR2oki+XilfUQvSZiGMlNU/TmoENsKBSUxIg0C06CXtroT32wB84F+7aaohN8B0FXEESbQKRdgwdBrO1V8iQ56FATWDoTLruKA02Ya16RmWr4uekZXygZh9SH0fE0DyJkd2jsdRzt4+Ex3uju6OEOXbIYSTcepIfqdMj6kGJwKtjSpnAduubLKY6/gOdKU+3UEdAPd548XGx9gDWoDufH+V87RlX/CC3GUByv0PxmTE/J1d8cBLehCJjD1clGm5uezUl0eRLFi76tQmgbG09tTkZC4xWjdicb6wNNQDjB8Zat9fyKTxz1qRKLmPnyyk37bqeuk6EETZgYkqnqB9eg2orpULR54Nh5tzRZWNfYoZEFuuZF6C8qd6K/KayIBpuAyR+7xw8AHKV1/UCJP0LchKb+RPFC6eEmNHalr+Z5orTyB68x/KfKSPj8mexlDO09TvyfBEumVelQ/HxXOt7rME6HB0vHWNUvISyh8OxBuxHnFdzSEvFEg/CJjFbOSdcQYEdkp91DRgnUpEPxQvTHJ4+dO4qM050dqeocoJRHiMo7QRbCzk7o/oQ0T5wSt8GE7bHAQXE3lB/pKVwMgxfV94YxIqepTZLr1nUoHthE+UQ1+gudqfOdU7/qHCDClRI52/qufOUgc07Ct7lonjzmjwbht5/aRXbiX7RTI8ZI0ns0JFCBDllDobcN9goakcD/XTRMxcq6e8AAAAAASUVORK5CYII=)
$\displaystyle - \kappa \left({\partial_{0} {\partial_{0} {{u}_{(0,0)}}}} + {\partial_{1} {\partial_{1} {{u}_{(0,0)}}}}\right) + \partial_t u_{C}$ $\displaystyle - \kappa \left({\partial_{0} {\partial_{0} {{u}_{(0,0)}}}} + {\partial_{1} {\partial_{1} {{u}_{(0,0)}}}}\right) + \partial_t u_{C}$
-κ⋅(D(D(u[0,0])) + D(D(u[0,0]))) + Transient(u_C) -κ⋅(D(D(u[0,0])) + D(D(u[0,0]))) + Transient(u_C)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Next, the PDE will be discretized. We use the `Discretization2ndOrder` class to apply finite differences discretization to the spatial components, and explicit euler discretization to the transient components. Next, the PDE will be discretized. We use the `Discretization2ndOrder` class to apply finite differences discretization to the spatial components, and explicit euler discretization to the transient components.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
discretize = ps.fd.Discretization2ndOrder(dx=dx, dt=dt) discretize = ps.fd.Discretization2ndOrder(dx=dx, dt=dt)
heat_pde_discretized = discretize(heat_pde) heat_pde_discretized = discretize(heat_pde)
heat_pde_discretized heat_pde_discretized
``` ```
%%%% Output: execute_result %%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAnkAAAAyCAYAAAAgCc0nAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAUBElEQVR4Ae2d67XUNhuFB9YpgEAHoQOSVBDoIJcKAh0ki3/5x/rSAUkFEDpIUgEJHYQOODkd8O3HxxKyxx5fxvZY9tZaHsm6v/uV5T26+c7Hjx8PNttE4Oeff/5Dkv0k+902JbRURsAIGAEjYAT2hYDe6fck8e+6vpX75pT0d0zyTsGTb5gUTwP4Q/avc0ihfD9Xvj+VeX8p+5p7+ZtQlqDYMgJ7RMB9wx61bpmXRkDP2SOVyXv+C7lbiZ5J3tKaWaA8KRzFv5cdSNikpSpfCN5L2U9CxnL/T+4fdT2R+8/gb9sIGIH9IOC+YT+6tqSXR0DP2zeqxf9kP2yrzd22APvniYCU/VQ1/1z2LASvRAVC9yxFqCyPfxMQzN0Zyf+4xH53svcR2Pj0QWkTcdw3JGp0u0/AaHEaoxZgengLuzeK9kZ263vXJK8HkLlEkaIZYaOT/XbmOj9W/v+qPNYFpIYRvHtlPVL/PbjBoo7HHuTuK6Px6YtU3vHcN1T153ZfxaPpzhg1odLTT+9bBnQeyWaA58iY5B1BkrUHGy1eSNnvZ5YCMsd08E1LOTy0NkbACOwPAfcN+9O5Jb48AhA9pm2P3r1Xl6+bazAFAiWLvy/7lynyO5WHymgbKWQh6EHhxeYL2TQ4/l2wXuAf3Vc2geieIeYfZLeRRQVv2xijbv0ao26M1hJDunLfMJEy3O67gTRGtxgJB6Zsn+vuaLmER/K621EuMVBuhUQtWXE1MAge08XpWsDn8od0MsJI/aKRP+TvG9mNBE/+hJHfbEb5s1Hk0ma1GK0EH/SzWozGNJ4V4Tqm+oPTSN6s+oYV6We17d4YDX4MeiU4E1fevU+VR+W96ZG8XtCvO1LZMBg1e3HBmjIqx7+JYiRRNh3727I+7MLliJXU4Nd43IrSQgCZDo7TznIHkvhBYQ90n5LJNN/o7kqj8F90vdRV2UQSM5jZoXKnxog28J2uZ8r7iz7VV7xWXBV2UXyov+owKUZ9MKnHUR02h2tdxpnvs+obttjuJVPrc96m+1Pt3hi1odbtPxeuyvdPXbxT0XUcUTfJ69ZJDjEgKRCsxlGxugBlI/tL/ry8+hoOXWwjZS+VCaQsNqzyPsSHeNQJKAu0636Hsm4cwxLzkhvS90F2IJDsZOUMwHiEi8IrZkCa3xWXtQydpJECFA9ZqXvd3MdD4U2E8Z38ozxJQjCbCiPIUKhXL72q7D64XhIfoJoMowT33k5htApce1f4zIiSl7bjvkGnBAiLS/ULaHGydt/zOa+0nJ7tftcYVQDreTMHrrWieaeiFzZAFnzA5+TVEMrttmw0/6jeHIgYCMNiYqhMiEKFlKWFKxziwXTtZ3IXjU42w8n/6jqqs8IgUTTSeNae3P/J72vZUT65+VTLQ9lxtE/30QxJo7jgR/5F/WImAxxKy3lFHF0zeE2k0pyNUahqWY+T5yYlcXvhqjwvig/1nRKjIP8Qe624DpFh6bjCLOu+QfXfRLuXHL2e86b20dXujVETat1+c+GqfPmThr75MEHxLrrbXR3HWDkCLLa8kUIjAVqqvmVDhWjFUSq5ITrpmgDC+FeaEigI0UF+TXV+LP+U4JEXDbdO5siPUZYjo/RD01BeGAU7ym8Bj7MwGlO/gRhdGh9EXByjneA6RszONGpfPOO59w3Zt/uBz3mnXhsiGKMGUCbwGoWr9M17kbRxRskkbwJtXDgLyAlKXdSoMUGwvpJdn+akc79OKgPhqpM5plmLOiv9U12QuINs4tJIU1OEpR6lmzLuN/jjNTQNawdbp35bypjS+1yMxtRlCEaXxgf5LoHRHnAdI+PJNHqOt9I3bKHdD3nOT+q1JdAYtQBzpvc5uP6usuNgy9WZFXHyCyJQdqY8xEyHLmZKMkZDYqEn06upYSQunbJkBC6SsbLOX8ov7ATm3/5NmQEv8pQglt6NFnm2dWCNCeTZlobyKftSZi6MxsjThNGl8UGONWG0JVzHyNKaRs8yz9FW+oYtt/um57xVrycCjNEJcM4IOgfX1yo3rB3/1STvDC2sIGmYYlx6JA9SSWfOmpu6qY/aMdL3mzp/dvx80MXLmp2frBvjCJNXuoK5JweNOzX1+xBGXPJqMkPT9CWWTWVN4XcuRmPqMASjS+ODfGdjpPb2SPmEXYanMGN3clvbOpWOsNxw7ZJnaPiW+oYttPsh7XGorom/W4xm7k9G46p6sXwLvTM7ZZJHK83YoEQUOvaFNEp0ldf6MeR6hmVji2v2kvAmP+S4l8Q5IFuZB6TyXRrWcF8Ej0hD3otimMpRyteER5PfEUZpXn3dAzG6KD7INAVGyoP2M+u0fG649m0vfeNJ/i31Ddm3+4Htsa+a03i7xWjm/uRcXP+WkvhTe7ibasvu4QhI0Xwzjq85/KerPnU5PMNhKZj2RJmbMMKPlzCNu27YFh5GLQ+Kh5sjYwpiJpvt4vURmpNpagVQJrt9zzH8c2r713xOvpW0krMNo0q8+s2ZGGWDD3KPxaiOWZ/7FeDap5rZxzmh05PP+Qr0s0i/gIInxqhvm9lK3zCmHfXFaEy8c3GN7wkfoTIG/oY0esA40oOz5N6EYLkZlWIbPkeMTDpSVObNVulf5X4WyszdlixM9xxNl8k/kDimfBktYIt4Qapk84+Fs734RFqKf2saxY1GaRrLjBFW5miqr/zoFGgHEGDwYF1kerbgaIyaylPeqzZT1dm4rkfNbTqVf+tzrjC3e6lwKEaKf7I/Ca1C8bLqO6l3W53lP7gdBRz62ipjEVxVDkupGHR6ctW3cluIJ8F5AbLrJCz6n0SsMl/yqq+NK8qT//UkBVUzobFgzh2Bus1lPb88aJCVyq5dYVy5T6urMP61cA4fO3ujOZUmRFIciPhB9qQkPOQ/k32EUVn/yTHKFB9gP8JojC5ywlV1naV/G4PbTGkadSq53e4/AT4JRl3tnuIUJ8e+k6pPghEZDTUL4hpm+B7dHVrJzOPTKIuGObEcrPMJa8fSrNv80zhj3YHk5UROOmXVQwBRZvo1yNeZhgiKzwuuTrL7pA0PfJ+4q4izMEbZ4YOSxmI0RsEj294cuM7Vv42BZfI0Y3W6Iv1Mjkk9Q2NUR+T4fmGMjivQ7TNF3xB4wcO9kbxueMfFaCMY+MfpwzRrNbRzyeb9Mr+gzDT7rN3ChpE8pmyHYMQI7c0QwRWfIW2+rpEdhqrz7BjljA/tYCRGQ5pQiDuo7eWOaxD6EvZIne5KP8aou2UugVF3LY5jTNU3KJ/wLrx/RTHy4GXKC4+1TmwiqExn6p5zj1jvFBLqdp9GGDDCBNOGGLA+jNEj1n280HVQOMSOFzCYEpeNGeD3VnZ6ftxfun8V/GQT/7ku9PAi+MvdZoiP2aROJD9r7oKMt5Ke+FXcSps9ETUNeq102eK3AEZZ44Oih2KUNo6+bpUxtO1lj2tfbOaIN1Sne9SPMepueQtg1F2J4xhT9w33io0XErb4ELNs1jRxptlnoWy5iwV8su8Ev1ztUj7+1aVkq7c4SgeBg7DxzdVi9Ec2C0/xj99mlfsgf7BklOgIN/lBXtg0wYYMDhSGJLJTFsLIYkkIDuvMWo3CIZqcM1cptzWBA4yAEdg0AuoT6HNG92+bBsfCGYGdIaD+AI7x/koOCMbbUn7WkNU3CeDXSjjKjuWd7FmmvJTvj7pGkbJSpkks1QFiBsGDgKWy4kb++ojQKdy+o1JKA8GjY+asuzAaQLo+hvrYGAEjYASMgBEwAkagjgBc7h7TtWwYCCQO8lFMOyaxGaWq+xXBSscoH+kj6ZGbESYMU5kPdN+686mIpZ9TaRT2i66XupgC7WWIr4jUu26KdWwteUHUmg6fDXkgFxsCAhkL/pTTtO6uzZ90EDmIHaOA1DXshJGztwmy1MnlyQxUJke92BgBI7BSBPSMHo3+p1VV+Bz928F9Q4qy3UZgfQh09Q0NNb7PSF5BEmRDShgdiiRGfqwpw+9o16LC8Ge6MRIjuSF96dlcj+X3h67W0ameaZj2LKaUlX+nUdxGQij/c6YzIMAVHJQfGIARZC2aNv8Y4ZaAvtY9I4OkZ30eRHn2dY8q5+QLRHWwMQJGYMUI6Bmeo387uG9YsdJdNSMwEoG7STrIWv0YEEgRD38Y6UuiF5sP+EeZGka7IhFSOtwQPYhMm+lMk+QDqVrcqHzK5aqQOd3HaddapYpRxLLeRVCZx0E2WJAXxJXpWg4z5luupIEkF0Z+XbIyFHvoEe82Q/8aASNgBIyAETACe0LgOiV5kI86mWMEriBtIhNPa4QC8pYSukBe4tRtieSNbNb9HRmlH5KmIIxHmSzrUZcNfArMJMs3pTzUKPpzQ5gsZMUcEcBb70+/ig/ZK6ZjP/keucAV0xXvNpZ/jYARMAJGwAgYgb0gADe4uUqkhcBEwiCi8Uj3X+oK07cP5VcQC9kQlkAy5CzMveCo2Yw4xXxrYUPSsDkE8tS0/q2W7bS3yK0LkhmI2kH3EDcwCkSXqetQN+QtCKH8kJEdbyEsEmf5100gkWAdcK/HCfesecS0YXgb6t9dI1C2v+clCKH9zr4sYNegW3gjkAEC7hsyUNL5VXyfkjw2SHB8CtOnEAgIB9OIrIX7UfYrXcHwsiimC4PHCRvCM5SINKW5UT7hJXWiuNmCmM4GH7B4oAs8IGwBn3TqGizZLFJMv8pOdweDRRpXt4X5Qb/fK+5XspvCb2N9+gUPGyPQhQDtM67hkpu29Y8uzsS0MQJGYL8IuG/Ytu7hGteR5KnzhzTETRSJ7E1+JK6TjPp9yIK4YYQq+AV7SJq+pDLkPak9BB/FRV4I4JFRWJs/I31htO8oXYNHwANCbLNBBNRWHkms33Tx54ZDMiNZGyAuyyyK9Z9lGv7E4cch3e8G5OOoRsAIrAQBnl9VxX3DSvSxtmqofcC7MDd3b+3Bv5CYkEmRWJniB2lrGm1rfJkMTEO+bWSxqEOPH+rHtQUTsGjCewvy7V4GPR8c68NoOs9afdNPX3wghn/3jex4WSOwpf4ta0XMXXn3DXMjnH3+gRd8iCN5Q0QqG1jIJE36QjePdRWkTvFwv5FdEBLZvKyey07PzjuZRvGDobx/w80YW+WG9XNjkq8qjWSBAFAnSIDNRhGQjnmGMKPartLX13ZC+tKzMYvM/ZM/AtL1qDaSv+T7lMB9wz713lPqwM/e3+2ZoCkaL4qQURGue9aecQByWKfGVC9rzYIhPlNFxdEsePZIE9IyzTlkOjOk27INeWZzjM12EaDd1482GiWtnjWmeHj2/MdgFIJOZARWhYD7hlWpY1WVCdzs3dUZ1WJtD6MC6agcpK1yn+avMEb4+NZqJHmEn0pThjMCSLwwRcmtze2IaRjpMR7bRAD9nj1Co2eHh55nlu8uM61nYwSMQN4IuG/IW39z1p4NnAVnGk3y9KLgIN9vdXE8SG/ypbhjGmYglHOCkmPerNPifD4+t+YXd44aTOosHQYixvPEDnfIHaNvLGmIBn3rJhyczcPMaDlpv9eFeas4cdS7zJdvLvPP/1DeY/d+bklnYwSMwGUQKJ9Z3oPuGy6jgtxKjTzr7jk1V8NjJO+Z7GKkrWdekMLehERxeZmxO9AvpGOAw4J6FGqTMQJq3+iQo00gY1wsfaBTx9RH8lgOwTedicP5keyye0w62TyLId1BfpA/jk3hSB921EIaiXety8YIGIGVI6Bn1n3DynW0puqVfT7vgWKz3lkkD8GUYev0bJPgil9fCN4ULfXj6Ij6Sy4N361buDD9DWEOIzi7xSJnwaVHHki+Ywy5S//M4GaDDToujNz86YkkTm7CWP4Qnqv7tXCIY3hJ4OZ6muapexsjYARWiICeU/cNK9TLyqsUBn2K2ZyrKSo75wtjzrynkH0FebxWHYpv6K6gLq7COAQgbUy5B6IWcuFhjdOupeffipcSQTZRRCKosMq5lrr/LGRm2wgYgewQcN+QncouXuGwIad4T5w9kndxcVwBRoAgCIG9G5H8EICkV0arpU/+wTPVWjkfT/7F8USJiKR9ldzbaQSMwHYQcN+wHV0uJQkzO/GrWSZ5S8E+Uzl66UMOYOyVEZyZinO2EyNQkjkIXYXM6b4YnS3121iqwiD2pI2jfeRX5tmYxp5GwAjkgUD5HLtvyENdq6il2kzYkBdnhSaZrl2FdPuuBEP6LKxnTdfNvqHIVvp0ChYhGHIvRu2kU/6Z4b7WxcgtGy8g92x8upE7TVs/bFxRbIyAEcgYgfT5Rgz3DRkrc+aqs0eCD1BEHuCRvJkRXyJ7KRTWjlKfL1Gey5gOgfJhhLAxNVsY+UHq2AUbdk8/kR8dPSN3XNe65x8+pC8a+RHGblsbI2AEMkdAzzN9uvuGzPW4VPXL/p/3SGUz7J2PHz8uVQeXMyMCUjDDtMzDc9h0ZPEzFumsJ0JA+oKwcQwKBO2BLtbYoUNGaPHjTEo2VxAPv+LzfrrnGBUIIf/s2TUL+YtTt7q3MQJGIGMEymfefUPGOlyq6morLPlhZqeydMskbykNLFCOlMuLHkJQYfILFO0ijIARMAJGwAgYgQsgoHc+sziQvKNBHk/XXkAhMxb5g/J+KoXHqb8Zy3LWRsAIGAEjYASMwOURYBaPL5AdzeKZ5F1eOZPVQApmcT6fwIrbpyfL3BkZASNgBIyAETACq0JA732W8LCcp3GpjkneqtR1fmWkaD51xbw8ircxAkbACBgBI2AENoiA3vPFZjzZlXV4qagmeSkaG3GXCudbpizKtzECRsAIGAEjYAQ2hIDe75zAwGDO16fEMsk7hU7GYWoAfO7qWdkQMpbEVTcCRsAIGAEjYAQCAnqvhxMZvpb7aB1eiIf9fyRZszVOHmoKAAAAAElFTkSuQmCC) ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAnkAAAAyCAYAAAAgCc0nAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAUBElEQVR4Ae2d67XUNhuFB9YpgEAHoQOSVBDoIJcKAh0ki3/5x/rSAUkFEDpIUgEJHYQOODkd8O3HxxKyxx5fxvZY9tZaHsm6v/uV5T26+c7Hjx8PNttE4Oeff/5Dkv0k+902JbRURsAIGAEjYAT2hYDe6fck8e+6vpX75pT0d0zyTsGTb5gUTwP4Q/avc0ihfD9Xvj+VeX8p+5p7+ZtQlqDYMgJ7RMB9wx61bpmXRkDP2SOVyXv+C7lbiZ5J3tKaWaA8KRzFv5cdSNikpSpfCN5L2U9CxnL/T+4fdT2R+8/gb9sIGIH9IOC+YT+6tqSXR0DP2zeqxf9kP2yrzd22APvniYCU/VQ1/1z2LASvRAVC9yxFqCyPfxMQzN0Zyf+4xH53svcR2Pj0QWkTcdw3JGp0u0/AaHEaoxZgengLuzeK9kZ263vXJK8HkLlEkaIZYaOT/XbmOj9W/v+qPNYFpIYRvHtlPVL/PbjBoo7HHuTuK6Px6YtU3vHcN1T153ZfxaPpzhg1odLTT+9bBnQeyWaA58iY5B1BkrUHGy1eSNnvZ5YCMsd08E1LOTy0NkbACOwPAfcN+9O5Jb48AhA9pm2P3r1Xl6+bazAFAiWLvy/7lynyO5WHymgbKWQh6EHhxeYL2TQ4/l2wXuAf3Vc2geieIeYfZLeRRQVv2xijbv0ao26M1hJDunLfMJEy3O67gTRGtxgJB6Zsn+vuaLmER/K621EuMVBuhUQtWXE1MAge08XpWsDn8od0MsJI/aKRP+TvG9mNBE/+hJHfbEb5s1Hk0ma1GK0EH/SzWozGNJ4V4Tqm+oPTSN6s+oYV6We17d4YDX4MeiU4E1fevU+VR+W96ZG8XtCvO1LZMBg1e3HBmjIqx7+JYiRRNh3727I+7MLliJXU4Nd43IrSQgCZDo7TznIHkvhBYQ90n5LJNN/o7kqj8F90vdRV2UQSM5jZoXKnxog28J2uZ8r7iz7VV7xWXBV2UXyov+owKUZ9MKnHUR02h2tdxpnvs+obttjuJVPrc96m+1Pt3hi1odbtPxeuyvdPXbxT0XUcUTfJ69ZJDjEgKRCsxlGxugBlI/tL/ry8+hoOXWwjZS+VCaQsNqzyPsSHeNQJKAu0636Hsm4cwxLzkhvS90F2IJDsZOUMwHiEi8IrZkCa3xWXtQydpJECFA9ZqXvd3MdD4U2E8Z38ozxJQjCbCiPIUKhXL72q7D64XhIfoJoMowT33k5htApce1f4zIiSl7bjvkGnBAiLS/ULaHGydt/zOa+0nJ7tftcYVQDreTMHrrWieaeiFzZAFnzA5+TVEMrttmw0/6jeHIgYCMNiYqhMiEKFlKWFKxziwXTtZ3IXjU42w8n/6jqqs8IgUTTSeNae3P/J72vZUT65+VTLQ9lxtE/30QxJo7jgR/5F/WImAxxKy3lFHF0zeE2k0pyNUahqWY+T5yYlcXvhqjwvig/1nRKjIP8Qe624DpFh6bjCLOu+QfXfRLuXHL2e86b20dXujVETat1+c+GqfPmThr75MEHxLrrbXR3HWDkCLLa8kUIjAVqqvmVDhWjFUSq5ITrpmgDC+FeaEigI0UF+TXV+LP+U4JEXDbdO5siPUZYjo/RD01BeGAU7ym8Bj7MwGlO/gRhdGh9EXByjneA6RszONGpfPOO59w3Zt/uBz3mnXhsiGKMGUCbwGoWr9M17kbRxRskkbwJtXDgLyAlKXdSoMUGwvpJdn+akc79OKgPhqpM5plmLOiv9U12QuINs4tJIU1OEpR6lmzLuN/jjNTQNawdbp35bypjS+1yMxtRlCEaXxgf5LoHRHnAdI+PJNHqOt9I3bKHdD3nOT+q1JdAYtQBzpvc5uP6usuNgy9WZFXHyCyJQdqY8xEyHLmZKMkZDYqEn06upYSQunbJkBC6SsbLOX8ov7ATm3/5NmQEv8pQglt6NFnm2dWCNCeTZlobyKftSZi6MxsjThNGl8UGONWG0JVzHyNKaRs8yz9FW+oYtt/um57xVrycCjNEJcM4IOgfX1yo3rB3/1STvDC2sIGmYYlx6JA9SSWfOmpu6qY/aMdL3mzp/dvx80MXLmp2frBvjCJNXuoK5JweNOzX1+xBGXPJqMkPT9CWWTWVN4XcuRmPqMASjS+ODfGdjpPb2SPmEXYanMGN3clvbOpWOsNxw7ZJnaPiW+oYttPsh7XGorom/W4xm7k9G46p6sXwLvTM7ZZJHK83YoEQUOvaFNEp0ldf6MeR6hmVji2v2kvAmP+S4l8Q5IFuZB6TyXRrWcF8Ej0hD3otimMpRyteER5PfEUZpXn3dAzG6KD7INAVGyoP2M+u0fG649m0vfeNJ/i31Ddm3+4Htsa+a03i7xWjm/uRcXP+WkvhTe7ibasvu4QhI0Xwzjq85/KerPnU5PMNhKZj2RJmbMMKPlzCNu27YFh5GLQ+Kh5sjYwpiJpvt4vURmpNpagVQJrt9zzH8c2r713xOvpW0krMNo0q8+s2ZGGWDD3KPxaiOWZ/7FeDap5rZxzmh05PP+Qr0s0i/gIInxqhvm9lK3zCmHfXFaEy8c3GN7wkfoTIG/oY0esA40oOz5N6EYLkZlWIbPkeMTDpSVObNVulf5X4WyszdlixM9xxNl8k/kDimfBktYIt4Qapk84+Fs734RFqKf2saxY1GaRrLjBFW5miqr/zoFGgHEGDwYF1kerbgaIyaylPeqzZT1dm4rkfNbTqVf+tzrjC3e6lwKEaKf7I/Ca1C8bLqO6l3W53lP7gdBRz62ipjEVxVDkupGHR6ctW3cluIJ8F5AbLrJCz6n0SsMl/yqq+NK8qT//UkBVUzobFgzh2Bus1lPb88aJCVyq5dYVy5T6urMP61cA4fO3ujOZUmRFIciPhB9qQkPOQ/k32EUVn/yTHKFB9gP8JojC5ywlV1naV/G4PbTGkadSq53e4/AT4JRl3tnuIUJ8e+k6pPghEZDTUL4hpm+B7dHVrJzOPTKIuGObEcrPMJa8fSrNv80zhj3YHk5UROOmXVQwBRZvo1yNeZhgiKzwuuTrL7pA0PfJ+4q4izMEbZ4YOSxmI0RsEj294cuM7Vv42BZfI0Y3W6Iv1Mjkk9Q2NUR+T4fmGMjivQ7TNF3xB4wcO9kbxueMfFaCMY+MfpwzRrNbRzyeb9Mr+gzDT7rN3ChpE8pmyHYMQI7c0QwRWfIW2+rpEdhqrz7BjljA/tYCRGQ5pQiDuo7eWOaxD6EvZIne5KP8aou2UugVF3LY5jTNU3KJ/wLrx/RTHy4GXKC4+1TmwiqExn6p5zj1jvFBLqdp9GGDDCBNOGGLA+jNEj1n280HVQOMSOFzCYEpeNGeD3VnZ6ftxfun8V/GQT/7ku9PAi+MvdZoiP2aROJD9r7oKMt5Ke+FXcSps9ETUNeq102eK3AEZZ44Oih2KUNo6+bpUxtO1lj2tfbOaIN1Sne9SPMepueQtg1F2J4xhT9w33io0XErb4ELNs1jRxptlnoWy5iwV8su8Ev1ztUj7+1aVkq7c4SgeBg7DxzdVi9Ec2C0/xj99mlfsgf7BklOgIN/lBXtg0wYYMDhSGJLJTFsLIYkkIDuvMWo3CIZqcM1cptzWBA4yAEdg0AuoT6HNG92+bBsfCGYGdIaD+AI7x/koOCMbbUn7WkNU3CeDXSjjKjuWd7FmmvJTvj7pGkbJSpkks1QFiBsGDgKWy4kb++ojQKdy+o1JKA8GjY+asuzAaQLo+hvrYGAEjYASMgBEwAkagjgBc7h7TtWwYCCQO8lFMOyaxGaWq+xXBSscoH+kj6ZGbESYMU5kPdN+686mIpZ9TaRT2i66XupgC7WWIr4jUu26KdWwteUHUmg6fDXkgFxsCAhkL/pTTtO6uzZ90EDmIHaOA1DXshJGztwmy1MnlyQxUJke92BgBI7BSBPSMHo3+p1VV+Bz928F9Q4qy3UZgfQh09Q0NNb7PSF5BEmRDShgdiiRGfqwpw+9o16LC8Ge6MRIjuSF96dlcj+X3h67W0ameaZj2LKaUlX+nUdxGQij/c6YzIMAVHJQfGIARZC2aNv8Y4ZaAvtY9I4OkZ30eRHn2dY8q5+QLRHWwMQJGYMUI6Bmeo387uG9YsdJdNSMwEoG7STrIWv0YEEgRD38Y6UuiF5sP+EeZGka7IhFSOtwQPYhMm+lMk+QDqVrcqHzK5aqQOd3HaddapYpRxLLeRVCZx0E2WJAXxJXpWg4z5luupIEkF0Z+XbIyFHvoEe82Q/8aASNgBIyAETACe0LgOiV5kI86mWMEriBtIhNPa4QC8pYSukBe4tRtieSNbNb9HRmlH5KmIIxHmSzrUZcNfArMJMs3pTzUKPpzQ5gsZMUcEcBb70+/ig/ZK6ZjP/keucAV0xXvNpZ/jYARMAJGwAgYgb0gADe4uUqkhcBEwiCi8Uj3X+oK07cP5VcQC9kQlkAy5CzMveCo2Yw4xXxrYUPSsDkE8tS0/q2W7bS3yK0LkhmI2kH3EDcwCkSXqetQN+QtCKH8kJEdbyEsEmf5100gkWAdcK/HCfesecS0YXgb6t9dI1C2v+clCKH9zr4sYNegW3gjkAEC7hsyUNL5VXyfkjw2SHB8CtOnEAgIB9OIrIX7UfYrXcHwsiimC4PHCRvCM5SINKW5UT7hJXWiuNmCmM4GH7B4oAs8IGwBn3TqGizZLFJMv8pOdweDRRpXt4X5Qb/fK+5XspvCb2N9+gUPGyPQhQDtM67hkpu29Y8uzsS0MQJGYL8IuG/Ytu7hGteR5KnzhzTETRSJ7E1+JK6TjPp9yIK4YYQq+AV7SJq+pDLkPak9BB/FRV4I4JFRWJs/I31htO8oXYNHwANCbLNBBNRWHkms33Tx54ZDMiNZGyAuyyyK9Z9lGv7E4cch3e8G5OOoRsAIrAQBnl9VxX3DSvSxtmqofcC7MDd3b+3Bv5CYkEmRWJniB2lrGm1rfJkMTEO+bWSxqEOPH+rHtQUTsGjCewvy7V4GPR8c68NoOs9afdNPX3wghn/3jex4WSOwpf4ta0XMXXn3DXMjnH3+gRd8iCN5Q0QqG1jIJE36QjePdRWkTvFwv5FdEBLZvKyey07PzjuZRvGDobx/w80YW+WG9XNjkq8qjWSBAFAnSIDNRhGQjnmGMKPartLX13ZC+tKzMYvM/ZM/AtL1qDaSv+T7lMB9wz713lPqwM/e3+2ZoCkaL4qQURGue9aecQByWKfGVC9rzYIhPlNFxdEsePZIE9IyzTlkOjOk27INeWZzjM12EaDd1482GiWtnjWmeHj2/MdgFIJOZARWhYD7hlWpY1WVCdzs3dUZ1WJtD6MC6agcpK1yn+avMEb4+NZqJHmEn0pThjMCSLwwRcmtze2IaRjpMR7bRAD9nj1Co2eHh55nlu8uM61nYwSMQN4IuG/IW39z1p4NnAVnGk3y9KLgIN9vdXE8SG/ypbhjGmYglHOCkmPerNPifD4+t+YXd44aTOosHQYixvPEDnfIHaNvLGmIBn3rJhyczcPMaDlpv9eFeas4cdS7zJdvLvPP/1DeY/d+bklnYwSMwGUQKJ9Z3oPuGy6jgtxKjTzr7jk1V8NjJO+Z7GKkrWdekMLehERxeZmxO9AvpGOAw4J6FGqTMQJq3+iQo00gY1wsfaBTx9RH8lgOwTedicP5keyye0w62TyLId1BfpA/jk3hSB921EIaiXety8YIGIGVI6Bn1n3DynW0puqVfT7vgWKz3lkkD8GUYev0bJPgil9fCN4ULfXj6Ij6Sy4N361buDD9DWEOIzi7xSJnwaVHHki+Ywy5S//M4GaDDToujNz86YkkTm7CWP4Qnqv7tXCIY3hJ4OZ6muapexsjYARWiICeU/cNK9TLyqsUBn2K2ZyrKSo75wtjzrynkH0FebxWHYpv6K6gLq7COAQgbUy5B6IWcuFhjdOupeffipcSQTZRRCKosMq5lrr/LGRm2wgYgewQcN+QncouXuGwIad4T5w9kndxcVwBRoAgCIG9G5H8EICkV0arpU/+wTPVWjkfT/7F8USJiKR9ldzbaQSMwHYQcN+wHV0uJQkzO/GrWSZ5S8E+Uzl66UMOYOyVEZyZinO2EyNQkjkIXYXM6b4YnS3121iqwiD2pI2jfeRX5tmYxp5GwAjkgUD5HLtvyENdq6il2kzYkBdnhSaZrl2FdPuuBEP6LKxnTdfNvqHIVvp0ChYhGHIvRu2kU/6Z4b7WxcgtGy8g92x8upE7TVs/bFxRbIyAEcgYgfT5Rgz3DRkrc+aqs0eCD1BEHuCRvJkRXyJ7KRTWjlKfL1Gey5gOgfJhhLAxNVsY+UHq2AUbdk8/kR8dPSN3XNe65x8+pC8a+RHGblsbI2AEMkdAzzN9uvuGzPW4VPXL/p/3SGUz7J2PHz8uVQeXMyMCUjDDtMzDc9h0ZPEzFumsJ0JA+oKwcQwKBO2BLtbYoUNGaPHjTEo2VxAPv+LzfrrnGBUIIf/s2TUL+YtTt7q3MQJGIGMEymfefUPGOlyq6morLPlhZqeydMskbykNLFCOlMuLHkJQYfILFO0ijIARMAJGwAgYgQsgoHc+sziQvKNBHk/XXkAhMxb5g/J+KoXHqb8Zy3LWRsAIGAEjYASMwOURYBaPL5AdzeKZ5F1eOZPVQApmcT6fwIrbpyfL3BkZASNgBIyAETACq0JA732W8LCcp3GpjkneqtR1fmWkaD51xbw8ircxAkbACBgBI2AENoiA3vPFZjzZlXV4qagmeSkaG3GXCudbpizKtzECRsAIGAEjYAQ2hIDe75zAwGDO16fEMsk7hU7GYWoAfO7qWdkQMpbEVTcCRsAIGAEjYAQCAnqvhxMZvpb7aB1eiIf9fyRZszVOHmoKAAAAAElFTkSuQmCC)
$\displaystyle {{u}_{(0,0)}} + dt \kappa \left(\frac{- 2 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(-1,0)}}}{dx^{2}} + \frac{- 2 {{u}_{(0,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}}}{dx^{2}}\right)$ $\displaystyle {{u}_{(0,0)}} + dt \kappa \left(\frac{- 2 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(-1,0)}}}{dx^{2}} + \frac{- 2 {{u}_{(0,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}}}{dx^{2}}\right)$
⎛-2⋅u_C + u_E + u_W -2⋅u_C + u_N + u_S⎞ ⎛-2⋅u_C + u_E + u_W -2⋅u_C + u_N + u_S⎞
u_C + dt⋅κ⋅⎜────────────────── + ──────────────────⎟ u_C + dt⋅κ⋅⎜────────────────── + ──────────────────⎟
⎜ 2 2 ⎟ ⎜ 2 2 ⎟
⎝ dx dx ⎠ ⎝ dx dx ⎠
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
It occurs to us that the right-hand summand can be simplified. It occurs to us that the right-hand summand can be simplified.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
heat_pde_discretized.simplify() heat_pde_discretized.simplify()
``` ```
%%%% Output: execute_result %%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAiUAAAAvCAYAAAA8abqkAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQJElEQVR4Ae2d6ZHVOBeGL1QH0AwZMBnAEMHXZMAMEQAZDMUv5h8FGQARsGQAEwFLBpABTGfQ3/uoLbft60WWr+/mV1VqW/vRo3OtY1l2X7u4uFjZdRP4559/TpX6tMhxqzg+VPx5d6nlpYjHI/X6h46fltd799gETMAETGCIgOaHV/KP+/Jd70t0WiDwQhCfFP5PxfyS/2o2VwTE5r5Cd3S0QXKFxWcmYAImYAJ1AsynH+tR9ZCNkjqPttAjQTyrJLzQ+S3F3a7ELfZUHMJKko691u9iAbnjJmACJmACgYDmiR86ea/j311IbJR0kbmKZ7L9chX0WYPAe4VfNeIcNAETMAETMIE1AjJIXivysY7c0K65k7WYI4pQp1nNeCPPXpB3Co++my8AVqlQB3snvlUjd3W+iT7myq624fqHjvdy60gtV/QTRR49hqltOJ8JmIAJmMBWCHAjy9zMloiaO+qVEk1g3+TvqMdYZL3PsWpUOgLFxBj2T3Rk2Xp0Xx+Vdir/XR7jYQ73RJVi9W7DsSLz2zYachsmYAImYAKzEmDeuM8c1WzlqI0SOqtOx/0gkzZhFhM7+0nY0HneBJkSRhZ53lLZqKPeosJmH4nHIGFz7hyOvkw29oYEU/86nz8OlT2W9Ll051j40A8z6h9N8+nnk5tqruPJiRlzKE8b/mqWPurHN0VnebTA45YsQ4I6VJaJnTdwwmOKIkw8m3bGOKxC/KZdVx+74ie3r77Hjb6z7rcp2mHsssdvcmf3o4K5dGc/ercZKcyon6P59PPJTTXXPHLMHTy+qa22H/1KiTrMakFzBUFRaU6TIgYJz794v/p2MUny2GKu1Yc0weq5uvpI/Id61suQ+sEPaYqj7knGXmLjDyRrTWkTyzmbCZiACZjA/hL4LtH+aIp3QkQxQbEU/7v81+YkoDDP8/f+g2GSEwOCRyysYPyUxxjhjv65fOkq/SXurvxDeco+kMd9Vp44mfNNEibw2rdJlL6TDZdDfVQ6xgKyITN9wpBi/OjTSx2j+1fhtzFOR/LzkTj04HmMj5lbjnA7b4nfWJRk4LFN55s9hcwHr7cbA9ZRkTl1gCmizaefD6lmNMwoJ8fCufL4hn2P+HIuuV6AfKpIJiz2BzCpl07xXPTZkFIWKhN1onjSmPxmcao7aT+B8jEZYzjED53Rn9iX5koJH3B5iVeez/LsAj5TmBUQJudYbqW4G/LXml55tu4kw2AfleeTPEtiYTLX+T3C8qVBonP6iLGGUqwU5pxne5RhSa3JS1FrjjpmWy2STOjUuY59j8j2Vm+hJdmTdJe8M7u95bQnjMxnWAH3ltGw6Os59kTvEGzJXOP8UbMfTjQ4TEhMzDj2IMSMIaKIC5NXjIhHlcVgYQm/nDh0Hid0VipuKsxE3+v6yigN4+GVfOfKhNKYIFkNwCApZVGYc97AOdcxOJ0jc5SRONJ4o4bVEtxv8tX0ELnrP2P6WMjKWLaOm+IxQFaqEwOGvjP5x0cklEtxcKqyTikzJg+v/3bqjtJ2rrd0RnJEXVnTd6UN6u4YIDl5Z+DEbw39YXx4s23Q7TOjGfh06kMXKMnQyVRpx6hDoxl1seuLXxrXPhZdaX2Mesp0jt9IfY3zMnNJ6Xh8U/3mBheb2qMOhbk7b8atis6Eu3ClB6c4JvyfCEaEjmfyH+U7JzqlpZThC3Dhc++hofU/QGIJKE6sMQeyx8cwMe6L8lUnUy6speGiNFYZJjnVwYoDbTddgK/0NgMLGfraHtNH2m3re5SH8cAQYWUsro7EtNQjF9KmARvKql7S/pXnmOpYzQlGlI7oRFjp6Sm8U71FrkLOIX0f0t1aF1XnpnVnk5wwBKNeJ43tATDaJJ+Ua1lzvFOYHpMOjWZUA5YYkN4timsillq2REbNMinjl6qvrfPHiQQ7p1UdudhwoSkndsWxrEJc23I+k2Rz4iDuf/LBqTx34hglfJa9agjELBwHyxT1YJRgeAR5qxXoHGOqJiN5FYf8tVdWFd9cPWgzxFQs36mNNqNjpXhWJWBRPkoZ0cqYPrb2vdIWY/1OntUlGLG/hPEZs2+IcQhGlo41p7pIS7qLrhVUQGWDzhXyNJPLcNEG+Xelt8iyCd0t+8SJ+rNR3VF9jAX1TuakOvjtYDyjx6lurxmpLxvjIyCDfW1CS2GqPFxH+65/tWqVd291KIdRrXOJgQVyTSRzlS2F0VXu8mxQx1Vvqr7G+aNmnFwvm7p8NYe7hvAjLeLDxacQvpI1nJ7ReIzUeZhMFG4aH9SH1brmRpahLS6sNac6ThWBrxkfCjOJr5Reyki46pQWL9Tlagr14av5dn1eyDOmj4GTypV9j33SMY4T1izK81oeA4IyWMHBKW5XDJDvrtpHvtIrLtz5FHH8MKJjdWmrekvDkiNyTNH3Vt2NHdjScRKnHBkPjNEkPiP7moPz4HVoC4wWyTWn0zllRo5fir7GOea8Ks9JJcBFtrmKwDJ/mNgkEBPWOx1Z9idvrSKFYwM6rTmsoGgR1RIUGFOGfS/IUxoQjcqakwN5Q38kL8YV58jC6gB3HfSLOwr6Uy3LxqPOvQzKv0tXlRM51vpY9KWMJ5Pi6D9lYbBmsCiu5pSfsYZPc4yr+aiva/yq+Uadq23aDTpXLaj4/4jXsfmIaxd6i2hdfUfHmvo+pLvUN7ebyilHvkNiNJXPmL7msDwGHZqb0VK55vQ7p8yY8UvR13id5JpZuuvlWWOFQxd/7kz/kGciw/2uuPPL0xU/4FpFRXzbgYa7OtOWn7i2MrRNuzVXyMQkVqYpjkkY+XmTBMfeFyZRJmT8L4WRqdYHxZEGzL1ykou+p/YR2eEXDJiinzwyiuOIwbI26VNILpTRkbGO5yGh5Q+c0I9tOcYL33Q1OSU3474rvUW2ZN1tdmTm8FyccsTeR0Zz8Wnraw4zrgHlNS6ngg2U2XdGOV08Zq45PHLKtOl4Clf0mUUB8pbupDzTmys6f6MMLI3zJgEKyLI+qwp/6/hWPrpTndQqaglX89aUOSaMLFMzICp1cMrdM7Ij5015ZGXyjbLHvS9MxuyZwfhYKT9vELDnhXReJ8ZY6VqJocguXWofkZGx5I2l8DhGx+oeFsYu8iBvdA91wofK7urYlh7zxSNGTnjUpTJNXYh5Jh9VN7LEi/GZwu8V/qxj7NMu9JZ+dfUZvk1979Nd6tqGm8opR8ZDYjSVz5i+5rA8Bh3KYqTf+m0Bqz6y7eLH9bz52+vKG+OPgetqZkaR1ZjxS+GKfbF2g1waJeoUDTaXxhGmLY6BPyUxOpShqIMJ5FuML47NcIgeWYZ6WxUuVfYiX20TmOIwQvbVECnwBQMqeXzgqoIYZWtOaV3xozioHh6lIBOrEmuKtdZwZoTaqI1Xs5pChjYdbYvbiN4ig9odo++dutvsz1zhqZxy5DokRlP5jOxrDs6D16FcRirH/NF63coB2Shz8Fzpz8yMArKR45fClblj7Qb4emht5J8CAI023XNFhFUIEpSP8w90pghzZ920eHvLUK5wtMdnaQ/ZMYnjj8XxBs9cF4uNM5LucXEbrbcIorK71t2t6U4Pp94xWQqjHj6917IOPr1MK4mbuP7tgw7NyaiCK/n0WLgmdzg1Y4e+9o5fpe5ertStvKyAMYfU3LWLi4taRGpAlfK2y9pymeKj0cEjID5bzwfNwkSsI0Lw/QpePS1XJ/rKKG9wytPaXkz3cfsE4njqeGP7ree12KVHiu/UW1pS+qJ0t42T4rjQPJbnZgMePEKrfqdlMYza+IjFSvGdeqS0NT6K62VKnTjlO7jrX5fMih/F6JLAuL9qY3FcxxEKOtXLSAzX9JU2+sYvyqA8vfqqdLYW1L5zFsuuMEpy/LNnz87kX2SWvT+mnNo5lf84pozz5o3rWG4al/fyo8ZzbBubzC9Zs/UWOcb2VfkPUnencFoCI/MZvr5sk1HONULyLe63mclp1PU9havyfJe/1SbP9dI6GXkiS4d9BCxpY20lO+XnLmvsHgQsa+7Q7PaPwEOJ9HT/xGqXKFdvqW1JupvLaSmMzKf991WN3TKjatOp5wc5r+RyTYVSzTfH71l1skrCixite0SzjRIEV6UYCjzC4flQquP11PA4J6WA8tIBPqLV2oGUOpxnPgLFWD7XMS7JztfYhmqWrDl6S+uL0t1MTothZD7DP8htMBqWYj2H5DroeSWT6zqI4ZiN/p4lN4sYtX8QuyZC2/LJ2DiWa8aWSc0/Z92pMjhf0lLtI43V2SGxmlu35q5/W6zn7MecdZvP8O/2GBjl9OEY9I5+71s/huRR+quhPNkbXdesG0eYgAmYgAmYgAmYwAQCkx7fTGjXRU3ABEzABEzABEygRuCallLy3gmuVeOACZiACZiACZiACUwj4Mc30/i5tAmYgAmYgAmYwIYI+PHNhkC6GhMwARMwARMwgWkEbJRM4+fSJmACJmACJmACGyJwsqF6XI0JmIAJ9BLQNwr4nlH80B7fK8DxLyeSv1t0WcR/TcAEjpWAjZJjHVn3ywT2j8CL4qNPQTKd8x9Cv8rzP7LsTMAETGDlxzdWAhMwgUECMiBuy3+V/68wJgbLtGR4pLL8m4no+AowX4zkH3/ZmYAJmMDKKyVWAhMwgUECMhy+KdMdHfmEAP8BNMfxef8vOQVdxgRMYBkEbJQsY5zdSxOYTKCyyjH2H2qGtlX+dUMIjJQfhcHTSHLQBExgiQT8+GaJo+4+m0AegXsqhhExeWOq6uCRzX35O3miuJQJmMAxErBRcoyj6j6ZwDwE2A+StUpSFUcGCW/esJ+Ex0GTDZxq3T43ARM4bAJ+fHPY42fpTWAWAhXD4Yca+CmPMcLqxnP50ikfr/nyb+Bxd+UfymN0PJDHfVaeD5enq1VR7xMdWXWJYY60Y2cCJrBwAjZKFq4A7r4JNAnIQGBF5L08KxnBWNAxbm5trpSUr/kqz98q80YeQwTDg1d+WREJRonCGCvEkRbfuGFfyRN5OxMwARPw2zfWARMwgSsCMhZY+cAgwXCorl5w/k1x5eMWnbNCgtERHWnsE2G1BPebfDWdb5JQP8fSqR4MEzsTMAETsFFiHTABE6gRwIg4laHQfFOG1ZPyMUxR4ovyVQ0XNq2WhovS/izyhYPCN6phn5uACZhAk4A3ujaJOGwCyybwl7pfe0QjY4LVDR69xEc4gZDi+XZJ1VH2bTXC5yZgAiYwhoCNkjG0nNcEjphAYXxggNSMD4UxNlZKrxkrxEWnNFZSKFuuplAfPubx0QRMwASGCHij6xAhp5vA8ghUH8nQe96UCasiMjLYM8L5L3n2nrDRFWOFfSHnOq+WfaqwN7EKjJ0JmEAaAa+UpHFyLhM4egIyIM7VSQwMHtUEpziMEN6U+XIZs7pXGB6sjOB/KcxqCEZK6RRH2ucywicmYAImkEDg2sUF/8rCzgRMwATCIxoMjPBar4435dkjgrHCBliMjE8yONjMSj7ivsuvFH4pjwHDqgpv12CslI9yFLYzARMwgUEC/weXMZOqdFgFVAAAAABJRU5ErkJggg==) ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAiUAAAAvCAYAAAA8abqkAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQJElEQVR4Ae2d6ZHVOBeGL1QH0AwZMBnAEMHXZMAMEQAZDMUv5h8FGQARsGQAEwFLBpABTGfQ3/uoLbft60WWr+/mV1VqW/vRo3OtY1l2X7u4uFjZdRP4559/TpX6tMhxqzg+VPx5d6nlpYjHI/X6h46fltd799gETMAETGCIgOaHV/KP+/Jd70t0WiDwQhCfFP5PxfyS/2o2VwTE5r5Cd3S0QXKFxWcmYAImYAJ1AsynH+tR9ZCNkjqPttAjQTyrJLzQ+S3F3a7ELfZUHMJKko691u9iAbnjJmACJmACgYDmiR86ea/j311IbJR0kbmKZ7L9chX0WYPAe4VfNeIcNAETMAETMIE1AjJIXivysY7c0K65k7WYI4pQp1nNeCPPXpB3Co++my8AVqlQB3snvlUjd3W+iT7myq624fqHjvdy60gtV/QTRR49hqltOJ8JmIAJmMBWCHAjy9zMloiaO+qVEk1g3+TvqMdYZL3PsWpUOgLFxBj2T3Rk2Xp0Xx+Vdir/XR7jYQ73RJVi9W7DsSLz2zYachsmYAImYAKzEmDeuM8c1WzlqI0SOqtOx/0gkzZhFhM7+0nY0HneBJkSRhZ53lLZqKPeosJmH4nHIGFz7hyOvkw29oYEU/86nz8OlT2W9Ll051j40A8z6h9N8+nnk5tqruPJiRlzKE8b/mqWPurHN0VnebTA45YsQ4I6VJaJnTdwwmOKIkw8m3bGOKxC/KZdVx+74ie3r77Hjb6z7rcp2mHsssdvcmf3o4K5dGc/ercZKcyon6P59PPJTTXXPHLMHTy+qa22H/1KiTrMakFzBUFRaU6TIgYJz794v/p2MUny2GKu1Yc0weq5uvpI/Id61suQ+sEPaYqj7knGXmLjDyRrTWkTyzmbCZiACZjA/hL4LtH+aIp3QkQxQbEU/7v81+YkoDDP8/f+g2GSEwOCRyysYPyUxxjhjv65fOkq/SXurvxDeco+kMd9Vp44mfNNEibw2rdJlL6TDZdDfVQ6xgKyITN9wpBi/OjTSx2j+1fhtzFOR/LzkTj04HmMj5lbjnA7b4nfWJRk4LFN55s9hcwHr7cbA9ZRkTl1gCmizaefD6lmNMwoJ8fCufL4hn2P+HIuuV6AfKpIJiz2BzCpl07xXPTZkFIWKhN1onjSmPxmcao7aT+B8jEZYzjED53Rn9iX5koJH3B5iVeez/LsAj5TmBUQJudYbqW4G/LXml55tu4kw2AfleeTPEtiYTLX+T3C8qVBonP6iLGGUqwU5pxne5RhSa3JS1FrjjpmWy2STOjUuY59j8j2Vm+hJdmTdJe8M7u95bQnjMxnWAH3ltGw6Os59kTvEGzJXOP8UbMfTjQ4TEhMzDj2IMSMIaKIC5NXjIhHlcVgYQm/nDh0Hid0VipuKsxE3+v6yigN4+GVfOfKhNKYIFkNwCApZVGYc97AOdcxOJ0jc5SRONJ4o4bVEtxv8tX0ELnrP2P6WMjKWLaOm+IxQFaqEwOGvjP5x0cklEtxcKqyTikzJg+v/3bqjtJ2rrd0RnJEXVnTd6UN6u4YIDl5Z+DEbw39YXx4s23Q7TOjGfh06kMXKMnQyVRpx6hDoxl1seuLXxrXPhZdaX2Mesp0jt9IfY3zMnNJ6Xh8U/3mBheb2qMOhbk7b8atis6Eu3ClB6c4JvyfCEaEjmfyH+U7JzqlpZThC3Dhc++hofU/QGIJKE6sMQeyx8cwMe6L8lUnUy6speGiNFYZJjnVwYoDbTddgK/0NgMLGfraHtNH2m3re5SH8cAQYWUsro7EtNQjF9KmARvKql7S/pXnmOpYzQlGlI7oRFjp6Sm8U71FrkLOIX0f0t1aF1XnpnVnk5wwBKNeJ43tATDaJJ+Ua1lzvFOYHpMOjWZUA5YYkN4timsillq2REbNMinjl6qvrfPHiQQ7p1UdudhwoSkndsWxrEJc23I+k2Rz4iDuf/LBqTx34hglfJa9agjELBwHyxT1YJRgeAR5qxXoHGOqJiN5FYf8tVdWFd9cPWgzxFQs36mNNqNjpXhWJWBRPkoZ0cqYPrb2vdIWY/1OntUlGLG/hPEZs2+IcQhGlo41p7pIS7qLrhVUQGWDzhXyNJPLcNEG+Xelt8iyCd0t+8SJ+rNR3VF9jAX1TuakOvjtYDyjx6lurxmpLxvjIyCDfW1CS2GqPFxH+65/tWqVd291KIdRrXOJgQVyTSRzlS2F0VXu8mxQx1Vvqr7G+aNmnFwvm7p8NYe7hvAjLeLDxacQvpI1nJ7ReIzUeZhMFG4aH9SH1brmRpahLS6sNac6ThWBrxkfCjOJr5Reyki46pQWL9Tlagr14av5dn1eyDOmj4GTypV9j33SMY4T1izK81oeA4IyWMHBKW5XDJDvrtpHvtIrLtz5FHH8MKJjdWmrekvDkiNyTNH3Vt2NHdjScRKnHBkPjNEkPiP7moPz4HVoC4wWyTWn0zllRo5fir7GOea8Ks9JJcBFtrmKwDJ/mNgkEBPWOx1Z9idvrSKFYwM6rTmsoGgR1RIUGFOGfS/IUxoQjcqakwN5Q38kL8YV58jC6gB3HfSLOwr6Uy3LxqPOvQzKv0tXlRM51vpY9KWMJ5Pi6D9lYbBmsCiu5pSfsYZPc4yr+aiva/yq+Uadq23aDTpXLaj4/4jXsfmIaxd6i2hdfUfHmvo+pLvUN7ebyilHvkNiNJXPmL7msDwGHZqb0VK55vQ7p8yY8UvR13id5JpZuuvlWWOFQxd/7kz/kGciw/2uuPPL0xU/4FpFRXzbgYa7OtOWn7i2MrRNuzVXyMQkVqYpjkkY+XmTBMfeFyZRJmT8L4WRqdYHxZEGzL1ykou+p/YR2eEXDJiinzwyiuOIwbI26VNILpTRkbGO5yGh5Q+c0I9tOcYL33Q1OSU3474rvUW2ZN1tdmTm8FyccsTeR0Zz8Wnraw4zrgHlNS6ngg2U2XdGOV08Zq45PHLKtOl4Clf0mUUB8pbupDzTmys6f6MMLI3zJgEKyLI+qwp/6/hWPrpTndQqaglX89aUOSaMLFMzICp1cMrdM7Ij5015ZGXyjbLHvS9MxuyZwfhYKT9vELDnhXReJ8ZY6VqJocguXWofkZGx5I2l8DhGx+oeFsYu8iBvdA91wofK7urYlh7zxSNGTnjUpTJNXYh5Jh9VN7LEi/GZwu8V/qxj7NMu9JZ+dfUZvk1979Nd6tqGm8opR8ZDYjSVz5i+5rA8Bh3KYqTf+m0Bqz6y7eLH9bz52+vKG+OPgetqZkaR1ZjxS+GKfbF2g1waJeoUDTaXxhGmLY6BPyUxOpShqIMJ5FuML47NcIgeWYZ6WxUuVfYiX20TmOIwQvbVECnwBQMqeXzgqoIYZWtOaV3xozioHh6lIBOrEmuKtdZwZoTaqI1Xs5pChjYdbYvbiN4ig9odo++dutvsz1zhqZxy5DokRlP5jOxrDs6D16FcRirH/NF63coB2Shz8Fzpz8yMArKR45fClblj7Qb4emht5J8CAI023XNFhFUIEpSP8w90pghzZ920eHvLUK5wtMdnaQ/ZMYnjj8XxBs9cF4uNM5LucXEbrbcIorK71t2t6U4Pp94xWQqjHj6917IOPr1MK4mbuP7tgw7NyaiCK/n0WLgmdzg1Y4e+9o5fpe5ertStvKyAMYfU3LWLi4taRGpAlfK2y9pymeKj0cEjID5bzwfNwkSsI0Lw/QpePS1XJ/rKKG9wytPaXkz3cfsE4njqeGP7ree12KVHiu/UW1pS+qJ0t42T4rjQPJbnZgMePEKrfqdlMYza+IjFSvGdeqS0NT6K62VKnTjlO7jrX5fMih/F6JLAuL9qY3FcxxEKOtXLSAzX9JU2+sYvyqA8vfqqdLYW1L5zFsuuMEpy/LNnz87kX2SWvT+mnNo5lf84pozz5o3rWG4al/fyo8ZzbBubzC9Zs/UWOcb2VfkPUnencFoCI/MZvr5sk1HONULyLe63mclp1PU9havyfJe/1SbP9dI6GXkiS4d9BCxpY20lO+XnLmvsHgQsa+7Q7PaPwEOJ9HT/xGqXKFdvqW1JupvLaSmMzKf991WN3TKjatOp5wc5r+RyTYVSzTfH71l1skrCixite0SzjRIEV6UYCjzC4flQquP11PA4J6WA8tIBPqLV2oGUOpxnPgLFWD7XMS7JztfYhmqWrDl6S+uL0t1MTothZD7DP8htMBqWYj2H5DroeSWT6zqI4ZiN/p4lN4sYtX8QuyZC2/LJ2DiWa8aWSc0/Z92pMjhf0lLtI43V2SGxmlu35q5/W6zn7MecdZvP8O/2GBjl9OEY9I5+71s/huRR+quhPNkbXdesG0eYgAmYgAmYgAmYwAQCkx7fTGjXRU3ABEzABEzABEygRuCallLy3gmuVeOACZiACZiACZiACUwj4Mc30/i5tAmYgAmYgAmYwIYI+PHNhkC6GhMwARMwARMwgWkEbJRM4+fSJmACJmACJmACGyJwsqF6XI0JmIAJ9BLQNwr4nlH80B7fK8DxLyeSv1t0WcR/TcAEjpWAjZJjHVn3ywT2j8CL4qNPQTKd8x9Cv8rzP7LsTMAETGDlxzdWAhMwgUECMiBuy3+V/68wJgbLtGR4pLL8m4no+AowX4zkH3/ZmYAJmMDKKyVWAhMwgUECMhy+KdMdHfmEAP8BNMfxef8vOQVdxgRMYBkEbJQsY5zdSxOYTKCyyjH2H2qGtlX+dUMIjJQfhcHTSHLQBExgiQT8+GaJo+4+m0AegXsqhhExeWOq6uCRzX35O3miuJQJmMAxErBRcoyj6j6ZwDwE2A+StUpSFUcGCW/esJ+Ex0GTDZxq3T43ARM4bAJ+fHPY42fpTWAWAhXD4Yca+CmPMcLqxnP50ikfr/nyb+Bxd+UfymN0PJDHfVaeD5enq1VR7xMdWXWJYY60Y2cCJrBwAjZKFq4A7r4JNAnIQGBF5L08KxnBWNAxbm5trpSUr/kqz98q80YeQwTDg1d+WREJRonCGCvEkRbfuGFfyRN5OxMwARPw2zfWARMwgSsCMhZY+cAgwXCorl5w/k1x5eMWnbNCgtERHWnsE2G1BPebfDWdb5JQP8fSqR4MEzsTMAETsFFiHTABE6gRwIg4laHQfFOG1ZPyMUxR4ovyVQ0XNq2WhovS/izyhYPCN6phn5uACZhAk4A3ujaJOGwCyybwl7pfe0QjY4LVDR69xEc4gZDi+XZJ1VH2bTXC5yZgAiYwhoCNkjG0nNcEjphAYXxggNSMD4UxNlZKrxkrxEWnNFZSKFuuplAfPubx0QRMwASGCHij6xAhp5vA8ghUH8nQe96UCasiMjLYM8L5L3n2nrDRFWOFfSHnOq+WfaqwN7EKjJ0JmEAaAa+UpHFyLhM4egIyIM7VSQwMHtUEpziMEN6U+XIZs7pXGB6sjOB/KcxqCEZK6RRH2ucywicmYAImkEDg2sUF/8rCzgRMwATCIxoMjPBar4435dkjgrHCBliMjE8yONjMSj7ivsuvFH4pjwHDqgpv12CslI9yFLYzARMwgUEC/weXMZOqdFgFVAAAAABJRU5ErkJggg==)
$\displaystyle \frac{{{u}_{(0,0)}} dx^{2} + dt \kappa \left(- 4 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}} + {{u}_{(-1,0)}}\right)}{dx^{2}}$ $\displaystyle \frac{{{u}_{(0,0)}} dx^{2} + dt \kappa \left(- 4 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}} + {{u}_{(-1,0)}}\right)}{dx^{2}}$
2 2
u_C⋅dx + dt⋅κ⋅(-4⋅u_C + u_E + u_N + u_S + u_W) u_C⋅dx + dt⋅κ⋅(-4⋅u_C + u_E + u_N + u_S + u_W)
─────────────────────────────────────────────── ───────────────────────────────────────────────
2 2
dx dx
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
While combining the two fractions on the right as desired, it also put everything above a common denominator. If we generated the kernel from this, we'd be redundantly multiplying $u_{(0,0)}$ by $dx^2$. Let's try something else. Instead of applying `simplify` to the entire equation, we could apply it only to the second summand. While combining the two fractions on the right as desired, it also put everything above a common denominator. If we generated the kernel from this, we'd be redundantly multiplying $u_{(0,0)}$ by $dx^2$. Let's try something else. Instead of applying `simplify` to the entire equation, we could apply it only to the second summand.
The outermost operation of `heat_pde_discretized` is a $+$, so `heat_pde_discretized` is an instance of `sp.Sum`. We take it apart by accessing its arguments, simplify the right hand summand, and put it back together again. The outermost operation of `heat_pde_discretized` is a $+$, so `heat_pde_discretized` is an instance of `sp.Sum`. We take it apart by accessing its arguments, simplify the right hand summand, and put it back together again.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify() heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify()
heat_pde_discretized heat_pde_discretized
``` ```
%%%% Output: execute_result %%%% Output: execute_result
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAgcAAAAvCAYAAAB3wWlJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQS0lEQVR4Ae2d6ZUctRqG23McwGBnABnYOII7ZOALEdjOwBz/gn8+kIEhAmxnAERgQwaQAb6Twdz3qSkJlWpTqbt6enl1jloq7Xpq+b7SUn3v5uZmY3M6BL7//vvn6s3fcn87nV65JyZgAiZgArsiIPnwRvbFVHkXU5GOOy4COtlP1eLHVgyO67y5tSZgAiawZwI/SE78OlXnPY8cTOE5njid6Eu19ne5j4+n1W6pCZiACZjAXRCQrGCU+VLuj0P1e+RgiMpxhr1Ts98cZ9PdahMwARMwgX0SkFLwk+p7IZcXy5653wtxQBUBAX6kjD/Lfi77VseT8zlVlYxkUl3U+aXcr0aS7Cy47ScX1N76t7PGuyATMAETMIGUAC+UyK3/poH4PXKQE6k8lrD8U5YhfbSwzlyOwhm6+UsWIb6G+VaFogXuwzBC8WAfFbkOEzABEzCBVQkgN54io/JarBzkRLY4FuCrNnu+U4BwFINPWxQ/lZW5o45CMpW4Nk79e1mb91TycY5l4W0zQsCMRsC0weYzzac21lyXkxOza+X6U/brPLeVg5zIdscM67ONEOCpGQtP01T5VRfTGZiPt846v2099Cvv2zoVHm6pl2oa1macgBmNsyHGfKb51Maaax05ZIenFerYFedihCAfNSAz4e/x5EZCd1tBQ9lDCkle1bbH36it+5q62Latzm8CJmACJlBG4C8l+zJPej8P8HEZAQlKpgl+kP1b9h9ZlALe4l/LbhSP0GbRHsKftI8Uxnz9B7np1hG2H/4SwuSS/pUsQ9evQ7j8Y+aJIq7HIncRrjYwnTC6E6JtM+39QvYPHXeUCB3T72dyV22n6jhoY07Tp8d8pvkQa0bzjGpSnDlXphVYF4eNz2hPK1RcSQKI4P9D9lv5sQh7FAVMM3KgsN9kGapphKr8X3EsGxUD+VEEUCg4ORsd42fuhzwM9QyNQii4YyhjrbUMtAnF5louStCYeaV4+sW6h8ChSatwlAYWvMSLrolofxROHHWsZlT+oayVOFhOB8LIfObvgoNlNN/0fooDue5o2DlzDfKj8xz2yEH/ep0M0cWMMOZNGKUgFZj42bGQC0HWGzTCX25umkUgyoMiwdcNEcLhrbt0W+ID5Uvbkdex7THbFtkNMWgUh0LzoY2kzeFCC+lH+6+8KA5MicT2yx+UC0ZjHup4tO5QAe5UPsX9KPtG9s62X6ruXXPiOuT64fwUffhK6UbZKu5OGa3AZ7SvXC9DRm0YZXrXfGjvITAa4jYXdm5c53gMxU8xGkpPmPKMXuOKW3I/B5mFLInGykFEUezhhDD8EoR4yMhowtC6grFw8iE4UQh44w6jBYQvMTzQcoHc5Fe5xP0ui1tqGN1olBm5CO/R6YS2QIR7UH4QVs20SlIZ/c/DNspDm5rRlJC2re8fuc3oitwr2V9lJxUlxdPOuXzvlI5PhpYqG/SbtuemuYFUzpCigXLYW9jTFrBLTigaoW1F51btOnRGu+RT0tfOeRWfEqandA0tZtQBVnhwblwLsXSSFTLK85Scv9LrdVB+3O/U6IMSAgjAznC/Ti4PaIZkOtsJx8KTSnjAv5VlJIL8rD/gLXrJ/Py10nc0Ph03RmURV/RW2WaJjvLSHpSg+FYfIxNPW8dGLn2BQ1SaFNaUobAOrzY7SlaueBD2nzaeMhlRQTn4XHaqHbP52rJQDugTXCaN0gwJ/43CGeGhPXF6aLKgNlLpmzrlbs1JZaCMoYjQllJz0Ix2yUdAZvuaQythqjRcjydxDdUwypmVHJ8h1xIsnTQljDoZbg9mr/EF12uQHx0lwcrBAPWxIMFG+GE7SoCO4/RAlhdBsOEkhXDKkGW0IAhOtLsQ/5P8/1NatMLw9lwkzEL5O3Rp3xO1B8UlNbxhIRwJ520vvInzxsxxKngb4aWwMLKQlsOoQBTA8gceuRJAedSZhzdlLcwHZ87J0AhPU94efrbiVNO+I2O0FZ+Ffa3BefTX0B4YnSXXmk7X5Fl4/kquV2QaJn12b6wc3EJZ+psLqjivrhOHQOStjjQxnAraOMIRlj3FgTSpUXqUBE5u56SlaeSnvHBys6j6Q9VNvdiOUTjKS1hsmcYh3HMlgP43ZbR94bPSQTHK+zTWB7TZoNmm9QX/knysjaBNd6kcbMsp9HuJe0yMtuWzpK9LGIa0p3ANrc0osFringLXJf3dJu2S81fCNTxfedZGcxF99swSQLApEcKOB1hjFIYywJstuwswzKMH5QHojV9hnFDeuIMAjYKTTJkJ+b9IysqSxENO/pfxaH0P/Ri6OEObmxao3TChXaG/9AV+GPh1LsQmdPgHhkP1Daf+N3QoH/XHc/dv0r361uJU04lDZLQWn6G+1jA75WtoV4zMtYbA9nmGzl/J9cozkZc20kbjkYOIotjDsOfPAvlS7kPZX2QR9MxFEpbOozPkzip5RgA2ctN56sssLUkwz2T54NATuWlZxA0ZhC9TD6tOP6h82hIEK1MC73T8QW7oE32FC3Nh7DTgIc96h8AFTsHQ986FOHCcpu0IjBDRunk5IZo68nylCkkoYw13W041bTomRtvyWdLXGpancA1VMdK9/UjAuL/nDDto8ntvLs8pcN2szCgwXHL+SrjynO6NEls5CLgLXZ18TgwKQm56Ye0NguLQM4obC2fIu3jYW+UwxE+beEvvneBexZUBqiOuDxgqom1Dj4HSDoXx4LhMy1H+sF4BBeTPNG7gOEYvzEfZSx9asa5deLblVNOGY2K0LZ+Ffa3BefTXUC0j5eO+HHxu1YDM8hw9V/qzMqMG2cLzV8IV2dF7Eb1oavPPsRN4qw6sddPunE17A3HR5ua1Apq1GEQoHf733AwcY+RnhCR/e5nNd5u7GfngU6HbGBQx7OpG/eRhPMRpsu5zYTTBZ/J6GOEzyTSJ5HycwjW0JqMEV7H3VLgWd7g04cj1Onn+krInuVK20jIihAzpmHs3NzedAB8cHwGdYE4u2yA/O5bWq63s+OgNPyo8CH6mJvgcMx+bisJY/qavCme7Zxxhmcun9BulGayTuEM1Q21WGDf8C1mUJ3gwtZN+5+FsGA3xEYuNwkevI8X1+ChskillYpTuJK6hti+LGDUAFv6cI9eFiLimJq89xfeu17nzF9qgvJPXq+KZ8u58bybk3aAcnIv97rvvrmSfn2J/1a93sk+PpW/tufihtr1L+6r0l7K/1tZ3V/m24XQOjMxn/vm9T0Y194nad3b3ZiWnRc/3Eq5K85fs50PtuYhawnl4GELBnqJ5pk69OpaOSWNlfQRTBGjNi4zy8Ma8dH0Fb0m8bR+VqeV0LozMZ/5y3jOj+Qb1U5zVvdnv/nzIGvezymTUgAXzcdo2bcm5KQdp30/KrxPM0PtruWGo8OD7p7YirJlaWKqwsSWU/hYZpeUm4GNTgzdBUSF3mEjtruF0NozMZ/7i3Aej+Vb0U6hd53hv9kHMh+z0fhZ3Xso6fwTYawLDCe3ww0u5b2R7w+4KY8j6cmjo4ZjC1Iensi+Pqc1L26r+PZe9WprvLtOrvateW2uXvy92a/ZjzbLNZ37o/xQY1fThFK47+n1o/Zhrj+KR9ZPP3WZBorQI9qKz8OuptAf2qseFbfKHoYd7Pc3iyALa/qGBhb35R9YDN9cETMAETMAE1idwIUHJSki+sodhO1z+0QTC8n3npG0MAld28bxxyD/nqmw+LGRjAiZgAiZgAiawJwL3VU/VX6XSPgluRhXIH+dy5Q9z3mxFe6hjvng2aabyKG7J/1I39SgPH3Rg0VpuHhCg+KGFaVN/t9uUo3ze95kT9bEJmIAJmMDJEYjfOZDgQ5iyJ/Iz+ZvFXnIZEeCDH4/l74we6JhFZExBxC/gyY+ywAr0ZtheLmUyXTH6gZ6SPG057MWcVTRU36hRfqZNPK0wSsgRJmACJmACJrDZpLsVEPLhE7aBDcJ0I6HaUQzaSEYI8k8uEha3mCkf/iu5U9MOs3mScpauam+bascETMAETMAETKCUQKocIMBzJYA3/kbYS0A/l02FM0I/VQTIT3ycYpAfcy3LuoaeUf4leRpFo1eIA0zABEzABEzABHZKgDUHwSDUmzl5AiS4Eej8IcNPHMvEv9xthTpCPzWp4pCGf9JBLDeNkH9JnpL/pc6K96EJmMCpE9DziOdI+AAYLxwYPq+dP6NuY/xrAiYwSyBVDpb8VSo3IEK/xKAYjCkBY/mH8nCjhxt/LJ/DTcAEzo8AW7HjImP5me78Q5b/5rAxAROoIBCVg1bLjosLk7KGwhD2uVaeH4ciSJtPNYS4JXlKlZFQtl0TMIEDJ6DnDiOUP8ui+L9NhfyCpjPlyRcwwzQn65gIeySbT5UuKNZJTeB8CVxUdh1hj9CPRjchYQj7obf7wRt0YR7KHVMyYjtmPLRvTCGZyepoEzCBXRNAeMs+Vrk8T9gtVWMYNfhYk9F5TMAEhgnEkYPh6OHQ9oYeUgJeK8eVbKMMKB3+93IboS6XB8AruemWxMk8Sh8M9W31P+qqN7xZhDLtmoAJ3DEB3Zc8JzBV96fyh3VRt6Xc/sFW+v2WEG7XBEygkMBFYbqhZNx8HQVBx3zfgA8fMQfIlw2ZkniWZCY9w33NFknC5Z/LE7Kzc+J9OLBrAiZwMgS4t/Nt1FWd0/OEaQqeL4xG2JiACVQSqBo5aOtiXo/hvHQUAGHfOU7bpThGFPjIUlQOiJ/K08Yz4kC6bacVKMbGBEzgsAgwclA1apB2Q88HXj54LvHRNk8fpnDsN4GFBKqVA918v8nyl498cbBYaCttzYMgKCILu+fkJmACh0SA54Xaw/3MM4NPrKMU8LbP9GI0SscLAV9cxTyRZQSSvN/IYj4oTRxJbMuNX2Ntjzdyi59Nt8X61wRMAAIX22DQjcfIwQu5zZt9YVkoE8VavdLygGAlsm/yQsBOZgKHSED3MC8GbDFEiGOZUkRRwOQjB0xN8r8qpOEbJ+xouCKfXJ43Id9GYSgNbF98Iz87FFA2SPdJ1sYETKCCwFbKAfXpRhydRhhqj9Lni4eGkqVhbG/KHxxpvP0mYAIHTkD3MAL9nSxKQaro42fHQnxhkJ8Xgij85SeOqcjw7HiQxaNwBMUDP5a1TbFMHduYgAksIFA9rZDWseZNuGbZaR/sNwETWJUAwp4/ZQsCPlSGUI/TA23gR6VLFQgWF0YFQnGdb6/o+LNQmF0TMIHdENh65GA3zXApJmACJ07ga/WvMwIooc5oAlMCne8bKDz/Lgp5f5G1MQET2BMBKwd7Au1qTOBcCbRKAIpARwnQMUJ/o/iO0kBYMIpjZIG8cXSB8rAhjV0TMIHdE9jJtMLum+USTcAETpBAOlVA9/i+QfhgGmsK8LOIkLUJLEhEaWDR87X8ad78Q2pKYmMCJrBLAh452CVNl2UCJtAjIMF+rUAEPVMIjVEYygC7Cj7ehmy+ahUARgqwn3TM6ADKQjQKI47dCzYmYAIrErh3c3OzYvEu2gRMwASaqQMEPdsREewPZVlDgNLAQkXC+G4Kiw5JR1jzqXQds50RRYJRBnYhoDTEKQYd25iACaxA4P8C3RkI5peT+QAAAABJRU5ErkJggg==) ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAgcAAAAvCAYAAAB3wWlJAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQS0lEQVR4Ae2d6ZUctRqG23McwGBnABnYOII7ZOALEdjOwBz/gn8+kIEhAmxnAERgQwaQAb6Twdz3qSkJlWpTqbt6enl1jloq7Xpq+b7SUn3v5uZmY3M6BL7//vvn6s3fcn87nV65JyZgAiZgArsiIPnwRvbFVHkXU5GOOy4COtlP1eLHVgyO67y5tSZgAiawZwI/SE78OlXnPY8cTOE5njid6Eu19ne5j4+n1W6pCZiACZjAXRCQrGCU+VLuj0P1e+RgiMpxhr1Ts98cZ9PdahMwARMwgX0SkFLwk+p7IZcXy5653wtxQBUBAX6kjD/Lfi77VseT8zlVlYxkUl3U+aXcr0aS7Cy47ScX1N76t7PGuyATMAETMIGUAC+UyK3/poH4PXKQE6k8lrD8U5YhfbSwzlyOwhm6+UsWIb6G+VaFogXuwzBC8WAfFbkOEzABEzCBVQkgN54io/JarBzkRLY4FuCrNnu+U4BwFINPWxQ/lZW5o45CMpW4Nk79e1mb91TycY5l4W0zQsCMRsC0weYzzac21lyXkxOza+X6U/brPLeVg5zIdscM67ONEOCpGQtP01T5VRfTGZiPt846v2099Cvv2zoVHm6pl2oa1macgBmNsyHGfKb51Maaax05ZIenFerYFedihCAfNSAz4e/x5EZCd1tBQ9lDCkle1bbH36it+5q62Latzm8CJmACJlBG4C8l+zJPej8P8HEZAQlKpgl+kP1b9h9ZlALe4l/LbhSP0GbRHsKftI8Uxnz9B7np1hG2H/4SwuSS/pUsQ9evQ7j8Y+aJIq7HIncRrjYwnTC6E6JtM+39QvYPHXeUCB3T72dyV22n6jhoY07Tp8d8pvkQa0bzjGpSnDlXphVYF4eNz2hPK1RcSQKI4P9D9lv5sQh7FAVMM3KgsN9kGapphKr8X3EsGxUD+VEEUCg4ORsd42fuhzwM9QyNQii4YyhjrbUMtAnF5louStCYeaV4+sW6h8ChSatwlAYWvMSLrolofxROHHWsZlT+oayVOFhOB8LIfObvgoNlNN/0fooDue5o2DlzDfKj8xz2yEH/ep0M0cWMMOZNGKUgFZj42bGQC0HWGzTCX25umkUgyoMiwdcNEcLhrbt0W+ID5Uvbkdex7THbFtkNMWgUh0LzoY2kzeFCC+lH+6+8KA5MicT2yx+UC0ZjHup4tO5QAe5UPsX9KPtG9s62X6ruXXPiOuT64fwUffhK6UbZKu5OGa3AZ7SvXC9DRm0YZXrXfGjvITAa4jYXdm5c53gMxU8xGkpPmPKMXuOKW3I/B5mFLInGykFEUezhhDD8EoR4yMhowtC6grFw8iE4UQh44w6jBYQvMTzQcoHc5Fe5xP0ui1tqGN1olBm5CO/R6YS2QIR7UH4QVs20SlIZ/c/DNspDm5rRlJC2re8fuc3oitwr2V9lJxUlxdPOuXzvlI5PhpYqG/SbtuemuYFUzpCigXLYW9jTFrBLTigaoW1F51btOnRGu+RT0tfOeRWfEqandA0tZtQBVnhwblwLsXSSFTLK85Scv9LrdVB+3O/U6IMSAgjAznC/Ti4PaIZkOtsJx8KTSnjAv5VlJIL8rD/gLXrJ/Py10nc0Ph03RmURV/RW2WaJjvLSHpSg+FYfIxNPW8dGLn2BQ1SaFNaUobAOrzY7SlaueBD2nzaeMhlRQTn4XHaqHbP52rJQDugTXCaN0gwJ/43CGeGhPXF6aLKgNlLpmzrlbs1JZaCMoYjQllJz0Ix2yUdAZvuaQythqjRcjydxDdUwypmVHJ8h1xIsnTQljDoZbg9mr/EF12uQHx0lwcrBAPWxIMFG+GE7SoCO4/RAlhdBsOEkhXDKkGW0IAhOtLsQ/5P8/1NatMLw9lwkzEL5O3Rp3xO1B8UlNbxhIRwJ520vvInzxsxxKngb4aWwMLKQlsOoQBTA8gceuRJAedSZhzdlLcwHZ87J0AhPU94efrbiVNO+I2O0FZ+Ffa3BefTX0B4YnSXXmk7X5Fl4/kquV2QaJn12b6wc3EJZ+psLqjivrhOHQOStjjQxnAraOMIRlj3FgTSpUXqUBE5u56SlaeSnvHBys6j6Q9VNvdiOUTjKS1hsmcYh3HMlgP43ZbR94bPSQTHK+zTWB7TZoNmm9QX/knysjaBNd6kcbMsp9HuJe0yMtuWzpK9LGIa0p3ANrc0osFringLXJf3dJu2S81fCNTxfedZGcxF99swSQLApEcKOB1hjFIYywJstuwswzKMH5QHojV9hnFDeuIMAjYKTTJkJ+b9IysqSxENO/pfxaH0P/Ri6OEObmxao3TChXaG/9AV+GPh1LsQmdPgHhkP1Daf+N3QoH/XHc/dv0r361uJU04lDZLQWn6G+1jA75WtoV4zMtYbA9nmGzl/J9cozkZc20kbjkYOIotjDsOfPAvlS7kPZX2QR9MxFEpbOozPkzip5RgA2ctN56sssLUkwz2T54NATuWlZxA0ZhC9TD6tOP6h82hIEK1MC73T8QW7oE32FC3Nh7DTgIc96h8AFTsHQ986FOHCcpu0IjBDRunk5IZo68nylCkkoYw13W041bTomRtvyWdLXGpancA1VMdK9/UjAuL/nDDto8ntvLs8pcN2szCgwXHL+SrjynO6NEls5CLgLXZ18TgwKQm56Ye0NguLQM4obC2fIu3jYW+UwxE+beEvvneBexZUBqiOuDxgqom1Dj4HSDoXx4LhMy1H+sF4BBeTPNG7gOEYvzEfZSx9asa5deLblVNOGY2K0LZ+Ffa3BefTXUC0j5eO+HHxu1YDM8hw9V/qzMqMG2cLzV8IV2dF7Eb1oavPPsRN4qw6sddPunE17A3HR5ua1Apq1GEQoHf733AwcY+RnhCR/e5nNd5u7GfngU6HbGBQx7OpG/eRhPMRpsu5zYTTBZ/J6GOEzyTSJ5HycwjW0JqMEV7H3VLgWd7g04cj1Onn+krInuVK20jIihAzpmHs3NzedAB8cHwGdYE4u2yA/O5bWq63s+OgNPyo8CH6mJvgcMx+bisJY/qavCme7Zxxhmcun9BulGayTuEM1Q21WGDf8C1mUJ3gwtZN+5+FsGA3xEYuNwkevI8X1+ChskillYpTuJK6hti+LGDUAFv6cI9eFiLimJq89xfeu17nzF9qgvJPXq+KZ8u58bybk3aAcnIv97rvvrmSfn2J/1a93sk+PpW/tufihtr1L+6r0l7K/1tZ3V/m24XQOjMxn/vm9T0Y194nad3b3ZiWnRc/3Eq5K85fs50PtuYhawnl4GELBnqJ5pk69OpaOSWNlfQRTBGjNi4zy8Ma8dH0Fb0m8bR+VqeV0LozMZ/5y3jOj+Qb1U5zVvdnv/nzIGvezymTUgAXzcdo2bcm5KQdp30/KrxPM0PtruWGo8OD7p7YirJlaWKqwsSWU/hYZpeUm4GNTgzdBUSF3mEjtruF0NozMZ/7i3Aej+Vb0U6hd53hv9kHMh+z0fhZ3Xso6fwTYawLDCe3ww0u5b2R7w+4KY8j6cmjo4ZjC1Iensi+Pqc1L26r+PZe9WprvLtOrvateW2uXvy92a/ZjzbLNZ37o/xQY1fThFK47+n1o/Zhrj+KR9ZPP3WZBorQI9qKz8OuptAf2qseFbfKHoYd7Pc3iyALa/qGBhb35R9YDN9cETMAETMAE1idwIUHJSki+sodhO1z+0QTC8n3npG0MAld28bxxyD/nqmw+LGRjAiZgAiZgAiawJwL3VU/VX6XSPgluRhXIH+dy5Q9z3mxFe6hjvng2aabyKG7J/1I39SgPH3Rg0VpuHhCg+KGFaVN/t9uUo3ze95kT9bEJmIAJmMDJEYjfOZDgQ5iyJ/Iz+ZvFXnIZEeCDH4/l74we6JhFZExBxC/gyY+ywAr0ZtheLmUyXTH6gZ6SPG057MWcVTRU36hRfqZNPK0wSsgRJmACJmACJrDZpLsVEPLhE7aBDcJ0I6HaUQzaSEYI8k8uEha3mCkf/iu5U9MOs3mScpauam+bascETMAETMAETKCUQKocIMBzJYA3/kbYS0A/l02FM0I/VQTIT3ycYpAfcy3LuoaeUf4leRpFo1eIA0zABEzABEzABHZKgDUHwSDUmzl5AiS4Eej8IcNPHMvEv9xthTpCPzWp4pCGf9JBLDeNkH9JnpL/pc6K96EJmMCpE9DziOdI+AAYLxwYPq+dP6NuY/xrAiYwSyBVDpb8VSo3IEK/xKAYjCkBY/mH8nCjhxt/LJ/DTcAEzo8AW7HjImP5me78Q5b/5rAxAROoIBCVg1bLjosLk7KGwhD2uVaeH4ciSJtPNYS4JXlKlZFQtl0TMIEDJ6DnDiOUP8ui+L9NhfyCpjPlyRcwwzQn65gIeySbT5UuKNZJTeB8CVxUdh1hj9CPRjchYQj7obf7wRt0YR7KHVMyYjtmPLRvTCGZyepoEzCBXRNAeMs+Vrk8T9gtVWMYNfhYk9F5TMAEhgnEkYPh6OHQ9oYeUgJeK8eVbKMMKB3+93IboS6XB8AruemWxMk8Sh8M9W31P+qqN7xZhDLtmoAJ3DEB3Zc8JzBV96fyh3VRt6Xc/sFW+v2WEG7XBEygkMBFYbqhZNx8HQVBx3zfgA8fMQfIlw2ZkniWZCY9w33NFknC5Z/LE7Kzc+J9OLBrAiZwMgS4t/Nt1FWd0/OEaQqeL4xG2JiACVQSqBo5aOtiXo/hvHQUAGHfOU7bpThGFPjIUlQOiJ/K08Yz4kC6bacVKMbGBEzgsAgwclA1apB2Q88HXj54LvHRNk8fpnDsN4GFBKqVA918v8nyl498cbBYaCttzYMgKCILu+fkJmACh0SA54Xaw/3MM4NPrKMU8LbP9GI0SscLAV9cxTyRZQSSvN/IYj4oTRxJbMuNX2Ntjzdyi59Nt8X61wRMAAIX22DQjcfIwQu5zZt9YVkoE8VavdLygGAlsm/yQsBOZgKHSED3MC8GbDFEiGOZUkRRwOQjB0xN8r8qpOEbJ+xouCKfXJ43Id9GYSgNbF98Iz87FFA2SPdJ1sYETKCCwFbKAfXpRhydRhhqj9Lni4eGkqVhbG/KHxxpvP0mYAIHTkD3MAL9nSxKQaro42fHQnxhkJ8Xgij85SeOqcjw7HiQxaNwBMUDP5a1TbFMHduYgAksIFA9rZDWseZNuGbZaR/sNwETWJUAwp4/ZQsCPlSGUI/TA23gR6VLFQgWF0YFQnGdb6/o+LNQmF0TMIHdENh65GA3zXApJmACJ07ga/WvMwIooc5oAlMCne8bKDz/Lgp5f5G1MQET2BMBKwd7Au1qTOBcCbRKAIpARwnQMUJ/o/iO0kBYMIpjZIG8cXSB8rAhjV0TMIHdE9jJtMLum+USTcAETpBAOlVA9/i+QfhgGmsK8LOIkLUJLEhEaWDR87X8ad78Q2pKYmMCJrBLAh452CVNl2UCJtAjIMF+rUAEPVMIjVEYygC7Cj7ehmy+ahUARgqwn3TM6ADKQjQKI47dCzYmYAIrErh3c3OzYvEu2gRMwASaqQMEPdsREewPZVlDgNLAQkXC+G4Kiw5JR1jzqXQds50RRYJRBnYhoDTEKQYd25iACaxA4P8C3RkI5peT+QAAAABJRU5ErkJggg==)
$\displaystyle {{u}_{(0,0)}} + \frac{dt \kappa \left(- 4 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}} + {{u}_{(-1,0)}}\right)}{dx^{2}}$ $\displaystyle {{u}_{(0,0)}} + \frac{dt \kappa \left(- 4 {{u}_{(0,0)}} + {{u}_{(1,0)}} + {{u}_{(0,1)}} + {{u}_{(0,-1)}} + {{u}_{(-1,0)}}\right)}{dx^{2}}$
dt⋅κ⋅(-4⋅u_C + u_E + u_N + u_S + u_W) dt⋅κ⋅(-4⋅u_C + u_E + u_N + u_S + u_W)
u_C + ───────────────────────────────────── u_C + ─────────────────────────────────────
2 2
dx dx
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
That looks a lot better! There is nothing left to simplify. The right-hand summand still contains a division by $dx^2$, though. Due to their inefficiency, floating-point divisions should be replaced by multiplication with their reciprocals. Before we can eliminate the division, we need to wrap the equation inside an `AssignmentCollection`. On this collection we can apply `add_subexpressions_for_divisions` to replace the division by a factor $\xi_1 = \frac{1}{dx^2}$ which in the kernel will be computed ahead of the loop. That looks a lot better! There is nothing left to simplify. The right-hand summand still contains a division by $dx^2$, though. Due to their inefficiency, floating-point divisions should be replaced by multiplication with their reciprocals. Before we can eliminate the division, we need to wrap the equation inside an `AssignmentCollection`. On this collection we can apply `add_subexpressions_for_divisions` to replace the division by a factor $\xi_1 = \frac{1}{dx^2}$ which in the kernel will be computed ahead of the loop.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
@ps.kernel @ps.kernel
def update(): def update():
u_tmp.center @= heat_pde_discretized u_tmp.center @= heat_pde_discretized
ac = ps.AssignmentCollection(update) ac = ps.AssignmentCollection(update)
ac = ps.simp.simplifications.add_subexpressions_for_divisions(ac) ac = ps.simp.simplifications.add_subexpressions_for_divisions(ac)
ac ac
``` ```
%%%% Output: execute_result %%%% Output: execute_result
AssignmentCollection: u_tmp_C, <- f(u_S, u_N, u_E, u_W, dt, u_C, kappa, dx) AssignmentCollection: u_tmp_C, <- f(u_S, u_N, u_E, u_W, dt, u_C, kappa, dx)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Our numeric solver's symbolic representation is now complete! Next, we use pystencils to generate and compile a C implementation of our kernel. The code is generated as shown below, compiled into a shared libary and then bound to `kernel_func`. All unbound sympy symbols (`dx`, `dt` and `kappa`) as well as the fields `u` and `u_tmp` are arguments to the generated kernel function. Our numeric solver's symbolic representation is now complete! Next, we use pystencils to generate and compile a C implementation of our kernel. The code is generated as shown below, compiled into a shared library and then bound to `kernel_func`. All unbound sympy symbols (`dx`, `dt` and `kappa`) as well as the fields `u` and `u_tmp` are arguments to the generated kernel function.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
kernel_ast = ps.create_kernel(update, cpu_openmp = 4) kernel_ast = ps.create_kernel(update, cpu_openmp = 4)
kernel_func = kernel_ast.compile() kernel_func = kernel_ast.compile()
ps.show_code(kernel_ast) ps.show_code(kernel_ast)
``` ```
%%%% Output: display_data %%%% Output: display_data
%%%% Output: display_data %%%% Output: display_data
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Prototype Simulation ### Prototype Simulation
We can set up and run a simple simulation wich the generated kernel right here. The first step is to set up the fields and simulation parameters. We can set up and run a simple simulation with the generated kernel right here. The first step is to set up the fields and simulation parameters.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
domain_size = 1.0 domain_size = 1.0
cells = 25 cells = 25
delta_x = domain_size / cells delta_x = domain_size / cells
delta_t = 0.0001 delta_t = 0.0001
kappa_v = 1.0 kappa_v = 1.0
u = np.zeros((cells, cells)) u = np.zeros((cells, cells))
u_tmp = np.zeros_like(u) u_tmp = np.zeros_like(u)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
We also need the Dirichlet and Neumann Boundaries. We also need the Dirichlet and Neumann Boundaries.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def f(x): def f(x):
return (1 + np.sin(2 * np.pi * x) * x**2) return (1 + np.sin(2 * np.pi * x) * x**2)
def init_domain(domain, domain_tmp): def init_domain(domain, domain_tmp):
domain.fill(0) domain.fill(0)
domain_tmp.fill(0) domain_tmp.fill(0)
domain[:,-1] = f( np.linspace(0, 1, domain.shape[0]) ) domain[:,-1] = f( np.linspace(0, 1, domain.shape[0]) )
domain_tmp[:,-1] = f( np.linspace(0, 1, domain_tmp.shape[0]) ) domain_tmp[:,-1] = f( np.linspace(0, 1, domain_tmp.shape[0]) )
return domain, domain_tmp return domain, domain_tmp