diff options
| author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2011-11-04 18:34:36 +0000 |
|---|---|---|
| committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2011-11-04 18:34:36 +0000 |
| commit | eb1d76b0841fa77d7158ec8392d5dcfbfba8949e (patch) | |
| tree | 456e08a4aac4f9fd3e478a456c9a6fe01f2dd79e | |
| parent | a2b91f57094795a9f825c5279f8962d52a7e3f8c (diff) | |
second iteration of constraint transport
* fix some indices
* move poison loop to proper spot
* compute conserved divergence in track_divB
* add Maple worksheet to check constraint transport indices
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@298 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
| -rw-r--r-- | doc/constraint_transport.mw | 268 | ||||
| -rw-r--r-- | schedule.ccl | 3 | ||||
| -rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 108 | ||||
| -rw-r--r-- | src/GRHydro_Source.F90 | 93 | ||||
| -rw-r--r-- | src/GRHydro_SourceM.F90 | 100 |
5 files changed, 429 insertions, 143 deletions
diff --git a/doc/constraint_transport.mw b/doc/constraint_transport.mw new file mode 100644 index 0000000..020e75c --- /dev/null +++ b/doc/constraint_transport.mw @@ -0,0 +1,268 @@ +<?xml version="1.0" encoding="UTF-8"?> +<Worksheet> +<Version major="15" minor="0"/> +<Label-Scheme value="2" prefix=""/> +<View-Properties presentation="false"><Zoom percentage="200"/></View-Properties> +<MapleNet-Properties elisiondigitsbefore="100" labelling="true" indentamount="4" elisiontermsthreshold="10000" ansi="false" errorbreak="1" useclientjvm="true" echo="1" imaginaryunit="I" labelwidth="20" contextmenusize="automatic" plotdriver="opengl" elisiondigitsafter="100" plotoutput="terminal" helpbrowser="standard" rtablesize="10" elisiontermsbefore="100" elisiondigitsthreshold="10000" typesetting="standard" plotdevice="inline" verboseproc="1" showassumed="1" errorcursor="false" longdelim="true" plotoptions="" quiet="false" elisiontermsafter="100" screenwidth="79" preplot="" prettyprint="3" displayprecision="-1" screenpixelheight="1080" warnlevel="3" screenheight="25" latexwidth="6.0" postplot="" prompt="> " ShowLabels="true"/> +<Styles><Font name="Ordered List 1" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Annotation Text" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Ordered List 2" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Ordered List 3" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Ordered List 4" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Ordered List 5" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Author" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Annotation Title" background="[255,255,255]" bold="true" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="18" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Warning" background="[255,255,255]" bold="false" executable="false" family="Courier New" foreground="[0,0,255]" italic="false" opaque="false" readonly="true" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Caption Reference" background="[255,255,255]" bold="true" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Maple Input Placeholder" background="[255,255,255]" bold="true" executable="true" family="Courier New" foreground="[200,0,200]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="true"/> +<Font name="Maple Plot" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Code" background="[255,255,255]" bold="false" executable="false" family="Courier New" foreground="[255,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Line Printed Output" background="[255,255,255]" bold="false" executable="false" family="Courier New" foreground="[0,0,255]" italic="false" opaque="false" readonly="true" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Text Output" background="[255,255,255]" bold="false" executable="false" family="Courier New" foreground="[0,0,255]" italic="false" opaque="false" readonly="true" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Diagnostic" background="[255,255,255]" bold="false" executable="false" family="Courier New" foreground="[40,120,40]" italic="false" opaque="false" readonly="true" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="2D Inert Output" background="[255,255,255]" bold="false" executable="true" family="Times New Roman" foreground="[144,144,144]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Normal" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Hyperlink" background="[255,255,255]" bold="false" executable="false" family="Times New Roma" foreground="[0,128,128]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="true" placeholder="false"/> +<Font name="Maple Output" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Dash Item" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="2D Math" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Maple Input" background="[0,0,0]" bold="true" executable="true" family="Monospaced" foreground="[255,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="2D Output" background="[0,0,0]" bold="false" executable="false" family="Lucida Bright" foreground="[0,0,255]" italic="false" opaque="false" readonly="true" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="2D Input" background="[255,255,255]" bold="false" executable="true" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="HyperlinkError" background="[255,255,255]" bold="false" executable="false" family="Courier New" foreground="[255,0,255]" italic="false" opaque="false" readonly="true" size="12" subscript="false" superscript="false" underline="true" placeholder="false"/> +<Font name="Header and Footer" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="10" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Error" background="[255,255,255]" bold="false" executable="false" family="Courier New" foreground="[255,0,255]" italic="false" opaque="false" readonly="true" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Title" background="[255,255,255]" bold="true" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="18" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Heading 1" background="[255,255,255]" bold="true" executable="false" family="Times New Roma" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="18" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Text" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Bullet Item" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Heading 4" background="[255,255,255]" bold="false" executable="false" family="Times New Roma" foreground="[0,0,0]" italic="true" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Equation Label" background="[255,255,255]" bold="true" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Heading 3" background="[255,255,255]" bold="true" executable="false" family="Times New Roma" foreground="[0,0,0]" italic="true" opaque="false" readonly="false" size="14" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="Heading 2" background="[255,255,255]" bold="true" executable="false" family="Times New Roma" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="16" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="HyperlinkWarning" background="[255,255,255]" bold="false" executable="false" family="Courier New" foreground="[0,0,255]" italic="false" opaque="false" readonly="true" size="12" subscript="false" superscript="false" underline="true" placeholder="false"/> +<Font name="Dictionary Hyperlink" background="[255,255,255]" bold="false" executable="false" family="Times New Roma" foreground="[147,0,15]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="true" placeholder="false"/> +<Font name="Caption Text" background="[255,255,255]" bold="true" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Font name="List Item" background="[255,255,255]" bold="false" executable="false" family="Times New Roman" foreground="[0,0,0]" italic="false" opaque="false" readonly="false" size="12" subscript="false" superscript="false" underline="false" placeholder="false"/> +<Layout name="Ordered List 1" alignment="left" bullet="numeric" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="3" spacebelow="3" linebreak="space" pagebreak-before="false" initial="-1" bulletsuffix=""/> +<Layout name="Ordered List 2" alignment="left" bullet="alphabetic" firstindent="0" leftmargin="36" rightmargin="0" linespacing="0.0" spaceabove="3" spacebelow="3" linebreak="space" pagebreak-before="false" initial="-1" bulletsuffix=""/> +<Layout name="Ordered List 3" alignment="left" bullet="roman" firstindent="0" leftmargin="72" rightmargin="0" linespacing="0.0" spaceabove="3" spacebelow="3" linebreak="space" pagebreak-before="false" initial="-1" bulletsuffix=""/> +<Layout name="Ordered List 4" alignment="left" bullet="ALPHABETIC" firstindent="0" leftmargin="108" rightmargin="0" linespacing="0.0" spaceabove="3" spacebelow="3" linebreak="space" pagebreak-before="false" initial="-1" bulletsuffix=""/> +<Layout name="Ordered List 5" alignment="left" bullet="ROMAN" firstindent="0" leftmargin="144" rightmargin="0" linespacing="0.0" spaceabove="3" spacebelow="3" linebreak="space" pagebreak-before="false" initial="-1" bulletsuffix=""/> +<Layout name="Author" alignment="centred" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="8" spacebelow="8" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Warning" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Annotation Title" alignment="centred" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="12" spacebelow="12" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Maple Plot" alignment="centred" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Line Printed Output" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="any" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Text Output" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="newline" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Diagnostic" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="any" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Normal" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Maple Output" alignment="centred" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.5" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Dash Item" alignment="left" bullet="dash" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="3" spacebelow="3" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="HyperlinkError" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Error" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Title" alignment="centred" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="12" spacebelow="12" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Heading 1" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="8" spacebelow="4" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Bullet Item" alignment="left" bullet="dot" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="3" spacebelow="3" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Heading 4" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Heading 3" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="Heading 2" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="8" spacebelow="2" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="HyperlinkWarning" alignment="left" bullet="none" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="0" spacebelow="0" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Layout name="List Item" alignment="left" bullet="indent" firstindent="0" leftmargin="0" rightmargin="0" linespacing="0.0" spaceabove="3" spacebelow="3" linebreak="space" pagebreak-before="false" initial="0" bulletsuffix=""/> +<Pencil-style name="Pencil 5" pen-color="[255,0,0]" pen-height="5.0" pen-width="5.0" pen-opacity="1.0"/> +<Pencil-style name="Pencil 4" pen-color="[0,0,255]" pen-height="3.0" pen-width="3.0" pen-opacity="1.0"/> +<Pencil-style name="Pencil 3" pen-color="[0,0,0]" pen-height="3.0" pen-width="3.0" pen-opacity="1.0"/> +<Pencil-style name="Pencil 2" pen-color="[0,0,255]" pen-height="1.0" pen-width="1.0" pen-opacity="1.0"/> +<Pencil-style name="Pencil 1" pen-color="[0,0,0]" pen-height="1.0" pen-width="1.0" pen-opacity="1.0"/> +<Highlighter-style name="Highlighter 2" pen-color="[255,204,0]" pen-height="14.0" pen-width="14.0" pen-opacity="0.8"/> +<Highlighter-style name="Highlighter 1" pen-color="[255,153,255]" pen-height="12.0" pen-width="8.0" pen-opacity="0.8"/> +<Highlighter-style name="Highlighter 4" pen-color="[0,255,255]" pen-height="32.0" pen-width="32.0" pen-opacity="0.8"/> +<Highlighter-style name="Highlighter 3" pen-color="[51,255,0]" pen-height="24.0" pen-width="24.0" pen-opacity="0.8"/> +<Highlighter-style name="Highlighter 5" pen-color="[255,255,0]" pen-height="48.0" pen-width="48.0" pen-opacity="0.8"/> +</Styles> +<Task-table> + <Task-category name="<default>"> + </Task-category> +</Task-table> +<Task> +</Task> +<Group labelreference="L34" drawlabel="true"> +<Input> +<Text-field style="Heading 1" layout="Heading 1"><Font executable="true">Constraint transport scheme in GRHydro</Font></Text-field> +<Text-field style="Text" layout="Heading 1">The current constraint transport scheme in GRHydro is based on WhiskyMHD (arXiv:gr-qc/0701109) and Bruno Giacomazzo's PhD thesis (http://www.sissa.it/ap/PhD/Theses/giacomazzo.ps.gz). It assumes that all magnetic field components live at the cell center and uses the Reconstruction step to move them to the cell interfaces.</Text-field> +<Text-field style="Text" layout="Heading 1"><Font encoding="UTF-8">This worksheet helps checking the expressions entered into the code by explicitely checking that div(\134dot B) = 0.</Font></Text-field> +</Input> +</Group> +<Group labelreference="L3"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">restart;</Text-field> +</Input> +</Group> +<Group labelreference="L4"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal"># we later expand taylor series in eps which helps counting orders in the stepsize +dx:=eps*dX;dy:=eps*dY;dz:=eps*dZ;</Text-field> +</Input> +</Group> +<Section collapsed="false" MultipleChoiceAnswerIndex="-1" MultipleChoiceRandomizeChoices="false" TrueFalseAnswerIndex="-1" EssayAnswerRows="5" EssayAnswerColumns="60"><Title> +<Text-field style="Heading 1" layout="Heading 1">Fluxes of the magnetic field components</Text-field></Title> +<Group labelreference="L5"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">BconsxfluxX := (i,j,k) -> 0;</Text-field> +</Input> +</Group> +<Group labelreference="L6"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">BconsyfluxX := (i,j,k) -> add(add(add(BconsyfluxX_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3);</Text-field> +</Input> +</Group> +<Group labelreference="L7"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">BconszfluxX := (i,j,k) -> add(add(add(BconszfluxX_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3);</Text-field> +</Input> +</Group> +<Group labelreference="L8"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal"># enforce F^I(B^j) = -F^j(B^i) if desired to check that RHS computed agrees to some order with the one used in divergence cleaning +BconsxfluxY := (i,j,k) -> add(add(add(BconsxfluxY_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); #-BconsyfluxX(i,j,k);</Text-field> +</Input> +</Group> +<Group labelreference="L9"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">BconsyfluxY := (i,j,k) -> 0;</Text-field> +</Input> +</Group> +<Group labelreference="L10"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">BconszfluxY := (i,j,k) -> add(add(add(BconszfluxY_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3);</Text-field> +</Input> +</Group> +<Group labelreference="L11"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">BconsxfluxZ := (i,j,k) -> add(add(add(BconsxfluxZ_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); #-BconszfluxX(i,j,k);</Text-field> +</Input> +</Group> +<Group labelreference="L12"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">BconsyfluxZ := (i,j,k) -> add(add(add(BconsyfluxZ_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); #-BconszfluxY(i,j,k);</Text-field> +</Input> +</Group> +<Group labelreference="L13"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">BconszfluxZ := (i,j,k) -> 0;</Text-field> +</Input> +</Group> +</Section> +<Section collapsed="false" MultipleChoiceAnswerIndex="-1" MultipleChoiceRandomizeChoices="false" TrueFalseAnswerIndex="-1" EssayAnswerRows="5" EssayAnswerColumns="60"><Title> +<Text-field style="Heading 1" layout="Heading 1">Lapse</Text-field></Title> +<Group labelreference="L14"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">alp := (i,j,k) -> add(add(add(alp_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3);</Text-field> +</Input> +</Group> +</Section> +<Section collapsed="false" MultipleChoiceAnswerIndex="-1" MultipleChoiceRandomizeChoices="false" TrueFalseAnswerIndex="-1" EssayAnswerRows="5" EssayAnswerColumns="60"><Title> +<Text-field style="Heading 1" layout="Heading 1"><Font encoding="UTF-8">Electric fields (eventually from E=v \134cross B=F(B))</Font></Text-field></Title> +<Group labelreference="L16"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">Ey := (i,j,k) -> 1/4*(alp(i+1/2,j,k)*BconszfluxX(i+1/2,j,k) + alp(i+1/2,j,k+1)*BconszfluxX(i+1/2,j,k+1)) - + 1/4*(alp(i,j,k+1/2)*BconsxfluxZ(i,j,k+1/2) + alp(i+1,j,k+1/2)*BconsxfluxZ(i+1,j,k+1/2));</Text-field> +</Input> +</Group> +<Group labelreference="L17"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">Ez := (i,j,k) -> 1/4*(alp(i,j+1/2,k)*BconsxfluxY(i,j+1/2,k) + alp(i+1,j+1/2,k)*BconsxfluxY(i+1,j+1/2,k)) - + 1/4*(alp(i+1/2,j,k)*BconsyfluxX(i+1/2,j,k) + alp(i+1/2,j+1,k)*BconsyfluxX(i+1/2,j+1,k));</Text-field> +</Input> +</Group> +<Group labelreference="L18"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">Ex := (i,j,k) -> 1/4*(alp(i,j,k+1/2)*BconsyfluxZ(i,j,k+1/2) + alp(i,j+1,k+1/2)*BconsyfluxZ(i,j+1,k+1/2)) - + 1/4*(alp(i,j+1/2,k)*BconszfluxY(i,j+1/2,k) + alp(i,j+1/2,k+1)*BconszfluxY(i,j+1/2,k+1));</Text-field> +</Input> +</Group> +</Section> +<Section collapsed="false" MultipleChoiceAnswerIndex="-1" MultipleChoiceRandomizeChoices="false" TrueFalseAnswerIndex="-1" EssayAnswerRows="5" EssayAnswerColumns="60"><Title> +<Text-field style="Heading 1" layout="Heading 1">RHS for B</Text-field></Title> +<Group labelreference="L19"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">Bconsrhsx := (i,j,k) -> -1/2*(Ey(i-1,j ,k-1)-Ey(i-1,j ,k ))/dz + -1/2*(Ez(i-1,j ,k )-Ez(i-1,j-1,k ))/dy + -1/2*(Ey(i ,j ,k-1)-Ey(i ,j ,k ))/dz + -1/2*(Ez(i ,j ,k )-Ez(i ,j-1,k ))/dy:</Text-field> +</Input> +</Group> +<Group labelreference="L20"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">Bconsrhsy := (i,j,k) -> -1/2 * (Ez(i-1,j-1,k )-Ez(i ,j-1,k )) / dx + -1/2 * (Ex(i ,j-1,k )-Ex(i ,j-1,k-1)) / dz + -1/2 * (Ez(i-1,j ,k )-Ez(i ,j ,k )) / dx + -1/2 * (Ex(i ,j ,k )-Ex(i ,j ,k-1)) / dz:</Text-field> +</Input> +</Group> +<Group labelreference="L21"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">Bconsrhsz := (i,j,k) -> - 1/2 * (Ex(i ,j-1,k-1)-Ex(i ,j ,k-1)) / dy + - 1/2 * (Ey(i ,j ,k-1)-Ey(i-1,j ,k-1)) / dx + - 1/2 * (Ex(i ,j-1,k )-Ex(i ,j ,k )) / dy + - 1/2 * (Ey(i ,j ,k )-Ey(i-1,j ,k )) / dx:</Text-field> +</Input> +</Group> +</Section> +<Section collapsed="false" MultipleChoiceAnswerIndex="-1" MultipleChoiceRandomizeChoices="false" TrueFalseAnswerIndex="-1" EssayAnswerRows="5" EssayAnswerColumns="60"><Title> +<Text-field style="Heading 1" layout="Heading 1">Check RHS at center of cell</Text-field></Title> +<Group labelreference="L36" drawlabel="true"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">Bconsrhsx(0,0,0):</Text-field> +</Input> +</Group> +<Group labelreference="L23"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">simplify(%):</Text-field> +</Input> +</Group> +<Group labelreference="L24"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">algsubs(eps=0,%);</Text-field> +</Input> +</Group> +<Group labelreference="L37" drawlabel="true"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal"></Text-field> +</Input> +</Group> +<Group labelreference="L40"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">(alp(0,-1/2,0)*BconsxfluxY(0,-1/2,0)-alp(0,1/2,0)*BconsxfluxY(0,1/2,0))/dy + +(alp(0,0,-1/2)*BconsxfluxZ(0,0,-1/2)-alp(0,0,1/2)*BconsxfluxZ(0,0,1/2))/dz:</Text-field> +</Input> +</Group> +<Group labelreference="L39"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">collect(%,eps):algsubs(eps^2=0,%);</Text-field> +</Input> +</Group> +</Section> +<Section collapsed="false" MultipleChoiceAnswerIndex="-1" MultipleChoiceRandomizeChoices="false" TrueFalseAnswerIndex="-1" EssayAnswerRows="5" EssayAnswerColumns="60"><Title> +<Text-field style="Heading 1" layout="Heading 1">Check constraint transport property for edge based divergence operator (see Bruno's thesis)</Text-field></Title> +<Group labelreference="L29"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">divB := (i,j,k) -> (Bconsrhsx(i+1,j,k)-Bconsrhsx(i,j,k)+Bconsrhsx(i+1,j+1,k)-Bconsrhsx(i,j+1,k)+Bconsrhsx(i+1,j,k+1)-Bconsrhsx(i,j,k+1)+Bconsrhsx(i+1,j+1,k+1)-Bconsrhsx(i,j+1,k+1))/4/dx + +(Bconsrhsy(i,j+1,k)-Bconsrhsy(i,j,k)+Bconsrhsy(i+1,j+1,k)-Bconsrhsy(i+1,j,k)+Bconsrhsy(i,j+1,k+1)-Bconsrhsy(i,j,k+1)+Bconsrhsy(i+1,j+1,k+1)-Bconsrhsy(i+1,j,k+1))/4/dy + +(Bconsrhsz(i,j,k+1)-Bconsrhsz(i,j,k)+Bconsrhsz(i+1,j,k+1)-Bconsrhsz(i+1,j,k)+Bconsrhsz(i,j+1,k+1)-Bconsrhsz(i,j+1,k)+Bconsrhsz(i+1,j+1,k+1)-Bconsrhsz(i+1,j+1,k))/4/dz;</Text-field> +</Input> +</Group> +<Group labelreference="L30"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal">simplify(divB(0,0,0));</Text-field> +</Input> +</Group> +</Section> +<Group labelreference="L41" drawlabel="true"> +<Input> +<Text-field prompt="> " style="Maple Input" layout="Normal"></Text-field> +</Input> +</Group> +<Text-field superscript="false" placeholder="false" executable="false" selection-placeholder="false" italic="false" size="12" bold="false" subscript="false" family="Times New Roman" opaque="false" underline="false" background="[255,255,255]" readonly="false" foreground="[0,0,0]" alignment="left" firstindent="0" spacebelow="0" leftmargin="0" linespacing="0.0" initial="0" linebreak="space" rightmargin="0" bulletsuffix="" spaceabove="0" bullet="none" pagebreak-before="false"></Text-field> +<Text-field superscript="false" placeholder="false" executable="false" selection-placeholder="false" italic="false" size="12" bold="false" subscript="false" family="Times New Roman" opaque="false" underline="false" background="[255,255,255]" readonly="false" foreground="[0,0,0]" alignment="left" firstindent="0" spacebelow="0" leftmargin="0" linespacing="0.0" initial="0" linebreak="space" rightmargin="0" bulletsuffix="" spaceabove="0" bullet="none" pagebreak-before="false"></Text-field> +</Worksheet>
\ No newline at end of file diff --git a/schedule.ccl b/schedule.ccl index ad870ce..311dcf4 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -107,6 +107,9 @@ if (evolve_tracer) STORAGE:GRHydro_tracer_rhs } +# FIXME: make temporary storage in FluxTerms or remove altogether, rename to be +# _not_ Evec which will likely collide with a future HydroBase electric field +# variable if (transport_constraints) { STORAGE:Evec diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90 index cc7b136..e53f4b8 100644 --- a/src/GRHydro_CalcUpdate.F90 +++ b/src/GRHydro_CalcUpdate.F90 @@ -39,7 +39,7 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) DECLARE_CCTK_FUNCTIONS CCTK_INT :: i,j,k,itracer - CCTK_REAL :: idx, alp_l, alp_r, Bvec_l, Bvec_r + CCTK_REAL :: idx, alp_l, alp_r, Bvec_l, Bvec_r, alp_tmp idx = 1.d0 / CCTK_DELTA_SPACE(flux_direction) @@ -73,39 +73,58 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * tauflux(i,j,k)) * idx if(evolve_mhd.ne.0) then - Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & - (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconsxflux(i,j,k)) * idx - Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + & - (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconsyflux(i,j,k)) * idx - Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + & - (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconszflux(i,j,k)) * idx + if(transport_constraints.ne.0) then + ! we have to first compute all components of v\crossB = E and + ! combine them in the last substep into Bconshs + ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k) + if(flux_direction.eq.1) then + alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) + Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) + alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) + Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) + elseif(flux_direction.eq.2) then + alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) + Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) + alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) + Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) + elseif(flux_direction.eq.3) then + alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) + Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) + alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) + Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) + end if + else + Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & + (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconsxflux(i,j,k)) * idx + Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + & + (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconsyflux(i,j,k)) * idx + Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + & + (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconszflux(i,j,k)) * idx + endif if(clean_divergence.ne.0) then psidcrhs(i,j,k) = psidcrhs(i,j,k) + & (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * psidcflux(i,j,k)) * idx endif - if(transport_constraints.ne.0) then - ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k) - if(flux_direction.eq.1) then - Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * alp_r*(Bconszflux(i,j,k) + Bconszflux(i,j ,k+1)) - Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsyflux(i,j+1,k )) - elseif(flux_direction.eq.2) then - Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * alp_r*(Bconszflux(i,j,k) + Bconszflux(i ,j,k+1)) - Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsxflux(i+1,j,k )) - elseif(flux_direction.eq.3) then - Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * alp_r*(Bconsyflux(i,j,k) + Bconsyflux(i ,j+1,k)) - Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsxflux(i+1,j ,k)) - end if - end if if(track_divB.ne.0) then - Bvec_l = 0.5d0 * (Bvec(i,j,k,flux_direction) + & - Bvec(i-xoffset,j-yoffset,k-zoffset,flux_direction)) - Bvec_r = 0.5d0 * (Bvec(i,j,k,flux_direction) + & - Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction)) - divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx + if(transport_constraints.ne.0) then + ! edge based divergence (see WhiskyMHD & Bruno's thesis, Eq. 7.27) + ! FIXME: this uses the Bcons before the current MoL substep update, should be in MoL_PseudoEvolution instead + divB(i,j,k) = divB(i,j,k) + & + 0.25d0*(Bcons(i+xoffset,j+yoffset ,k+zoffset,flux_direction)-Bcons(i ,j ,k ,flux_direction)+ & + Bcons(i+1 ,j+1-zoffset,k+zoffset,flux_direction)-Bcons(i+1-xoffset,j+xoffset ,k ,flux_direction)+ & + Bcons(i+xoffset,j+1-xoffset,k+1 ,flux_direction)-Bcons(i ,j+zoffset ,k+1-zoffset,flux_direction)+ & + Bcons(i+1 ,j+1 ,k+1 ,flux_direction)-Bcons(i+1-xoffset,j+1-yoffset,k+1-zoffset,flux_direction))*idx + else + Bvec_l = 0.5d0 * (Bvec(i,j,k,flux_direction) + & + Bvec(i-xoffset,j-yoffset,k-zoffset,flux_direction)) + Bvec_r = 0.5d0 * (Bvec(i,j,k,flux_direction) + & + Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction)) + divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx + endif endif endif @@ -272,32 +291,27 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) end if if (transport_constraints.ne.0 .and. flux_direction.eq.1) then ! HACK: x direction is last - ! RH: I think one could wrap all of this into a single do loop and remove the - ! Evec storage + ! FIXME: I think one could wrap all of this into a single do loop and remove the + ! Evec storage do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil - Bconsrhs(i,j,k,1) = & !Bconsrhs(i,j,k,1) & - - 0.5d0 * (Evec(i-1,j-1,k-1,2)-Evec(i-1,j-1,k ,2)) / CCTK_DELTA_SPACE(3) & - - 0.5d0 * (Evec(i-1,j ,k-1,3)-Evec(i-1,j-1,k-1,3)) / CCTK_DELTA_SPACE(2) & - - 0.5d0 * (Evec(i ,j-1,k-1,2)-Evec(i ,j-1,k ,2)) / CCTK_DELTA_SPACE(3) & - - 0.5d0 * (Evec(i ,j ,k-1,3)-Evec(i ,j-1,k-1,3)) / CCTK_DELTA_SPACE(2) - Bconsrhs(i,j,k,2) = & !Bconsrhs(i,j,k,2) & - - 0.5d0 * (Evec(i-1,j-1,k-1,3)-Evec(i-1,j ,k-1,3)) / CCTK_DELTA_SPACE(1) & - - 0.5d0 * (Evec(i ,j-1,k-1,1)-Evec(i-1,j-1,k-1,1)) / CCTK_DELTA_SPACE(3) & - - 0.5d0 * (Evec(i-1,j-1,k ,3)-Evec(i-1,j ,k ,3)) / CCTK_DELTA_SPACE(1) & - - 0.5d0 * (Evec(i ,j-1,k ,1)-Evec(i-1,j-1,k ,1)) / CCTK_DELTA_SPACE(3) - Bconsrhs(i,j,k,3) = & !Bconsrhs(i,j,k,3) & - - 0.5d0 * (Evec(i-1,j-1,k-1,1)-Evec(i ,j-1,k-1,1)) / CCTK_DELTA_SPACE(1) & - - 0.5d0 * (Evec(i-1,j-1,k ,2)-Evec(i-1,j-1,k-1,2)) / CCTK_DELTA_SPACE(2) & - - 0.5d0 * (Evec(i-1,j ,k-1,1)-Evec(i ,j ,k-1,1)) / CCTK_DELTA_SPACE(1) & - - 0.5d0 * (Evec(i-1,j ,k ,2)-Evec(i-1,j ,k-1,2)) / CCTK_DELTA_SPACE(2) + Bconsrhs(i,j,k,1) = - 0.5d0 * (Evec(i-1,j ,k-1,2)-Evec(i-1,j ,k ,2)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i-1,j ,k ,3)-Evec(i-1,j-1,k ,3)) / CCTK_DELTA_SPACE(2) & + - 0.5d0 * (Evec(i ,j ,k-1,2)-Evec(i ,j ,k ,2)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i ,j ,k ,3)-Evec(i ,j-1,k ,3)) / CCTK_DELTA_SPACE(2) + Bconsrhs(i,j,k,2) = - 0.5d0 * (Evec(i-1,j-1,k ,3)-Evec(i ,j-1,k ,3)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i ,j-1,k ,1)-Evec(i ,j-1,k-1,1)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i-1,j ,k ,3)-Evec(i ,j ,k ,3)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i ,j ,k ,1)-Evec(i ,j ,k-1,1)) / CCTK_DELTA_SPACE(3) + Bconsrhs(i,j,k,3) = - 0.5d0 * (Evec(i ,j-1,k-1,1)-Evec(i ,j ,k-1,1)) / CCTK_DELTA_SPACE(2) & + - 0.5d0 * (Evec(i ,j ,k-1,2)-Evec(i-1,j ,k-1,2)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i ,j-1,k ,1)-Evec(i ,j ,k ,1)) / CCTK_DELTA_SPACE(2) & + - 0.5d0 * (Evec(i ,j ,k ,2)-Evec(i-1,j ,k ,2)) / CCTK_DELTA_SPACE(1) enddo enddo enddo end if - - return diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90 index 7b04b8c..c72ea5e 100644 --- a/src/GRHydro_Source.F90 +++ b/src/GRHydro_Source.F90 @@ -454,99 +454,6 @@ subroutine SourceTerms(CCTK_ARGUMENTS) deallocate(force_spatial_second_order) -#if(0) /* poison edges of domain */ - if(last_iteration_seen .ne. cctk_iteration .or. reflevel .ne. grhydro_reflevel) then - last_iteration_seen = cctk_iteration - reflevel = grhydro_reflevel - mol_substep = 0 - else - mol_substep = mol_substep + 1 - end if - do k = 1, GRHydro_stencil*mol_substep - do j = 1, cctk_lsh(2) - do i = 1, cctk_lsh(1) - dens(i,j,k) = -1d100 - Scon(i,j,k,1) = -1d100 - Scon(i,j,k,2) = -1d100 - Scon(i,j,k,3) = -1d100 - tau(i,j,k) = -1d100 - Bcons(i,j,k,1) = -1d100 - Bcons(i,j,k,2) = -1d100 - Bcons(i,j,k,3) = -1d100 - end do - end do - end do - do k = cctk_lsh(3)-GRHydro_stencil*mol_substep+1, cctk_lsh(3) - do j = 1, cctk_lsh(2) - do i = 1, cctk_lsh(1) - dens(i,j,k) = -1d100 - Scon(i,j,k,1) = -1d100 - Scon(i,j,k,2) = -1d100 - Scon(i,j,k,3) = -1d100 - tau(i,j,k) = -1d100 - Bcons(i,j,k,1) = -1d100 - Bcons(i,j,k,2) = -1d100 - Bcons(i,j,k,3) = -1d100 - end do - end do - end do - do i = 1, GRHydro_stencil*mol_substep - do k = 1, cctk_lsh(3) - do j = 1, cctk_lsh(2) - dens(i,j,k) = -1d100 - Scon(i,j,k,1) = -1d100 - Scon(i,j,k,2) = -1d100 - Scon(i,j,k,3) = -1d100 - tau(i,j,k) = -1d100 - Bcons(i,j,k,1) = -1d100 - Bcons(i,j,k,2) = -1d100 - Bcons(i,j,k,3) = -1d100 - end do - end do - end do - do i = cctk_lsh(1)-GRHydro_stencil*mol_substep+1, cctk_lsh(1) - do k = 1, cctk_lsh(3) - do j = 1, cctk_lsh(2) - dens(i,j,k) = -1d100 - Scon(i,j,k,1) = -1d100 - Scon(i,j,k,2) = -1d100 - Scon(i,j,k,3) = -1d100 - tau(i,j,k) = -1d100 - Bcons(i,j,k,1) = -1d100 - Bcons(i,j,k,2) = -1d100 - Bcons(i,j,k,3) = -1d100 - end do - end do - end do - do j = 1, GRHydro_stencil*mol_substep - do i = 1, cctk_lsh(1) - do k = 1, cctk_lsh(3) - dens(i,j,k) = -1d100 - Scon(i,j,k,1) = -1d100 - Scon(i,j,k,2) = -1d100 - Scon(i,j,k,3) = -1d100 - tau(i,j,k) = -1d100 - Bcons(i,j,k,1) = -1d100 - Bcons(i,j,k,2) = -1d100 - Bcons(i,j,k,3) = -1d100 - end do - end do - end do - do j = cctk_lsh(2)-GRHydro_stencil*mol_substep+1, cctk_lsh(2) - do i = 1, cctk_lsh(1) - do k = 1, cctk_lsh(3) - dens(i,j,k) = -1d100 - Scon(i,j,k,1) = -1d100 - Scon(i,j,k,2) = -1d100 - Scon(i,j,k,3) = -1d100 - tau(i,j,k) = -1d100 - Bcons(i,j,k,1) = -1d100 - Bcons(i,j,k,2) = -1d100 - Bcons(i,j,k,3) = -1d100 - end do - end do - end do -#endif end subroutine SourceTerms diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90 index 1456500..cdf9ef0 100644 --- a/src/GRHydro_SourceM.F90 +++ b/src/GRHydro_SourceM.F90 @@ -113,6 +113,9 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) Bconsrhs=0.d0 if (evolve_tracer .ne. 0) then + Bconsrhsx(i,j,k) = 0.d0 + Bconsrhsy(i,j,k) = 0.d0 + Bconsrhsz(i,j,k) = 0.d0 cons_tracerrhs = 0.d0 end if @@ -130,9 +133,6 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) if (transport_constraints .ne. 0) then Evec = 0.d0 - Bconsrhsx(i,j,k) = 0.d0 - Bconsrhsy(i,j,k) = 0.d0 - Bconsrhsz(i,j,k) = 0.d0 endif @@ -495,6 +495,100 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) deallocate(force_spatial_second_order) +#if(0) /* poison edges of domain */ + if(last_iteration_seen .ne. cctk_iteration .or. reflevel .ne. grhydro_reflevel) then + last_iteration_seen = cctk_iteration + reflevel = grhydro_reflevel + mol_substep = 0 + else + mol_substep = mol_substep + 1 + end if + do k = 1, GRHydro_stencil*mol_substep + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do k = cctk_lsh(3)-GRHydro_stencil*mol_substep+1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do i = 1, GRHydro_stencil*mol_substep + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do i = cctk_lsh(1)-GRHydro_stencil*mol_substep+1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do j = 1, GRHydro_stencil*mol_substep + do i = 1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do j = cctk_lsh(2)-GRHydro_stencil*mol_substep+1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do +#endif + end subroutine SourceTermsM |