<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoAcetate, li.MsoAcetate, div.MsoAcetate
{mso-style-priority:99;
mso-style-link:"Balloon Text Char";
margin:0cm;
margin-bottom:.0001pt;
font-size:8.0pt;
font-family:"Tahoma","sans-serif";
mso-fareast-language:EN-US;}
p.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0cm;
margin-right:0cm;
margin-bottom:0cm;
margin-left:36.0pt;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
span.EmailStyle18
{mso-style-type:personal;
font-family:"Calibri","sans-serif";
color:windowtext;}
span.EmailStyle19
{mso-style-type:personal;
font-family:"Calibri","sans-serif";
color:#1F497D;}
span.EmailStyle20
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
span.BalloonTextChar
{mso-style-name:"Balloon Text Char";
mso-style-priority:99;
mso-style-link:"Balloon Text";
font-family:"Tahoma","sans-serif";
mso-fareast-language:EN-US;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:72.0pt 72.0pt 72.0pt 72.0pt;}
div.WordSection1
{page:WordSection1;}
/* List Definitions */
@list l0
{mso-list-id:1847477534;
mso-list-type:hybrid;
mso-list-template-ids:1034564988 -355715036 134807577 134807579 134807567 134807577 134807579 134807567 134807577 134807579;}
@list l0:level1
{mso-level-text:"\(%1\)";
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;}
@list l0:level2
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;}
@list l0:level3
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level4
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;}
@list l0:level5
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;}
@list l0:level6
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level7
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;}
@list l0:level8
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-18.0pt;}
@list l0:level9
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
ol
{margin-bottom:0cm;}
ul
{margin-bottom:0cm;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-GB" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal"><span style="color:#1F497D">Dear All,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D">In case anyone ever has the same sort issue in the future I thought I would submit the solution to the mailing list for archiving.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D">This problem came about from an ambiguity in my given attempt resulting from a failure to explicitly distinguish between gridDataStruct == CENTER or FACEX,FACEY,FACEZ<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D">If you look at Flash.h in the object directory, you see that the cell-centred variables are numbered independently of the face variables (e.g., CURZ_VAR=1, DENS_VAR=2,.. etc, and MAG_FACE_VAR=1,MAGI_FACE_VAR=2).
Failing to distinguish between the two allowed the “if ivar = MAGI_FACE_VAR” statement to be evaluated as true for the density also.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D">If I had looked more carefully at the examples for the outflow/diode/reflecting boundary conditions, I would have noticed that this distinction is actually made when selecting which sign = +-1 is applied. Schoolboy
error!<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D">Cheers,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D">Jonathan<o:p></o:p></span></p>
<p class="MsoNormal"><a name="_MailEndCompose"><span style="color:#1F497D"><o:p> </o:p></span></a></p>
<div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0cm 0cm 0cm">
<p class="MsoNormal"><b><span lang="EN-US" style="font-size:10.0pt;font-family:"Tahoma","sans-serif";mso-fareast-language:EN-GB">From:</span></b><span lang="EN-US" style="font-size:10.0pt;font-family:"Tahoma","sans-serif";mso-fareast-language:EN-GB"> flash-users-bounces@flash.uchicago.edu
[mailto:flash-users-bounces@flash.uchicago.edu] <b>On Behalf Of </b>Jonathan Thurgood<br>
<b>Sent:</b> 05 May 2016 13:41<br>
<b>To:</b> flash-users@flash.uchicago.edu<br>
<b>Subject:</b> [FLASH-USERS] User Boundary conditions for USM - WARNING after gc filling: min. unk(DENS_VAR)=-1<o:p></o:p></span></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><span style="color:black">Dear All, <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">I spotted some of the sillier errors I sent you all in my message about problems at a 2D null point yesterday (apologies) and got further with the problem.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">I now suffer from an error message indicating that the boundary conditions (BC) I attempted somehow forces a –ve density:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">(e.g. WARNING after gc filling: min. unk(DENS_VAR)=-1.000000000000000000000 PE=0 block=1 type=1 )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">I’ve been able to reproduce this error in a simplified problem of a 2D homogeneous plasma equilibrium (dens = 1 eint=1e-12 vel=0) with a uniform field of Bx = -1, By = 0.
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">I attempt to set the B values in the boundaries to their correct value of B=[-1,0] “manually” and get the above error, indicating that I have somehow assigned the value intended for the field to the density.
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">The relevant bits of my Grid_bcApplyToRegionSpecialized.F90 file are attached below, and commented where appropriate.
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal" style="margin-bottom:12.0pt"><span style="color:black">I have
<o:p></o:p></span></p>
<p class="MsoListParagraph" style="margin-bottom:12.0pt;text-indent:-18.0pt;mso-list:l0 level1 lfo2">
<![if !supportLists]><span style="color:black"><span style="mso-list:Ignore">(1)<span style="font:7.0pt "Times New Roman"">
</span></span></span><![endif]><span style="color:black">set only the magx_var and magy_var (cell-centered) manually whilst changing nothing else – this works<o:p></o:p></span></p>
<p class="MsoListParagraph" style="margin-bottom:12.0pt;text-indent:-18.0pt;mso-list:l0 level1 lfo2">
<![if !supportLists]><span style="color:black"><span style="mso-list:Ignore">(2)<span style="font:7.0pt "Times New Roman"">
</span></span></span><![endif]><span style="color:black">set only the MAG_FACE_VAR values manually – this produces a running code, but introduces an artefact which propagates in from the boundary of the order 1e-7 <o:p></o:p></span></p>
<p class="MsoListParagraph" style="text-indent:-18.0pt;mso-list:l0 level1 lfo2"><![if !supportLists]><span style="color:black"><span style="mso-list:Ignore">(3)<span style="font:7.0pt "Times New Roman"">
</span></span></span><![endif]><span style="color:black">set only the MAGI_FACE_VAR values manually - this line seems to be responsible for forcing a –ve density In particular. If this value is changed, say from the “real” value of -1 to -2, the error will
report that gc filling set DENS_VAR=given negative value there.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">Could anyone let me know why (2) introduces the artefact (and if possible how to avoid), and why (3) is apparently setting values for variables other than MAGI_FACE_VAR?
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">More generally, any extra advice about setting user BC (particularly for USM) would be appreciated also. As far as I know, for a custom BC in USM I need only to set velocities (velx_var, vely_var, velz_var), cell-centred
fields (magx_var, magy_var, magz_var), dens, eint (or other depending on EOS) and face-fields (MAG_FACE_VAR and MAGI_FACE_VAR). <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">Thank you,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">Jonathan<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">!!!!!!!!! Bits of Code Below here !!!!!!!!<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"> do ivar = 1,varCount<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if(mask(ivar)) then<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> call gr_bcMapBcType(bcTypeActual,bcType,ivar,gridDataStruct,axis,face,idest)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if(face==LOW) then<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> select case (bcTypeActual)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> case(USER_DEFINED)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> k = 2*guard+1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if(isFace)k=k+1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> do i = 1,guard<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> regionData(i,1:je,1:ke,ivar) = regionData(k-i,1:je,1:ke,ivar) !all variables are uniform and in equilibrium - this should set
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> !correct BC values for all variables<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"> !as test, then going to overwrite this to set some variables (namely, fields) manually<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> !to try and pinpoint the error<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> !the idea is eventually to set the field in the boundary as f(x,y)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">!1. the following two specifications for the cell-centre fields work as expected - no artefacts introduced<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">!<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">#ifdef MAGX_VAR <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if( ivar==MAGX_VAR ) regionData(i,1:je,1:ke,ivar) = -1.0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">#endif<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">#ifdef MAGY_VAR<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if( ivar==MAGY_VAR ) regionData(i,1:je,1:ke,ivar) = 0.0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">#endif<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">!2. The following specifies bx on the i-drn face as -1 and by on the j-drn face as 0.
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! the code will run with these lines un-commented but with
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! a small propagating artefact is launched from the lower corner<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! error in magx in domain is of order 1e-7<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">#ifdef MAG_FACE_VAR<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if( (ivar==MAG_FACE_VAR).AND.(axis==IAXIS) ) regionData(i,1:je,1:ke,ivar) = -1.0
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if( (ivar==MAG_FACE_VAR).AND.(axis==JAXIS) ) regionData(i,1:je,1:ke,ivar) = 0.0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">#endif<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">!3. Same thing with MAGI<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! This produces a crash as follows<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> ! Zero or imaginary sound speed has obtained! Please try other (more diffusive) slope limiter, flux, order, cfl, etc.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! WARNING after gc filling: min. unk(DENS_VAR)=-1.000000000000000000000 PE=0 block=1 type=1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! WARNING after gc filling: min. unk(DENS_VAR)=0.000000000000000000000 PE=4 block=1 type=1
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! (etc for other blocks)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">!<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">! Looks as if this has actually set dens_var = -1 or =0 rather than magi_face_var<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">!<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black">#ifdef MAGI_FACE_VAR<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if( (ivar==MAGI_FACE_VAR).AND.(axis==IAXIS) ) regionData(i,1:je,1:ke,ivar) = -1.0
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if( (ivar==MAGI_FACE_VAR).AND.(axis==JAXIS) ) regionData(i,1:je,1:ke,ivar) = 0.0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">#endif<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> end do<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> end select<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> else !face==HIGH<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"> ! (I only tested the above in the lower-facing boundaries)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> select case (bcTypeActual)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> case(USER_DEFINED)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> k = 2*guard+1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> if(isFace)k=k+1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> do i = 1,guard<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> regionData(k-i,1:je,1:ke,ivar)= regionData(i,1:je,1:ke,ivar)
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> end do<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> end select<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> endif<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"> endif <o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"> enddo !end ivar loop<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:black"> return<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black">end subroutine Grid_bcApplyToRegionSpecialized<o:p></o:p></span></p>
<p class="MsoNormal"><span style="color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="color:#1F497D"><o:p> </o:p></span></p>
<div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0cm 0cm 0cm">
<p class="MsoNormal"><b><span lang="EN-US" style="font-size:10.0pt;font-family:"Tahoma","sans-serif";mso-fareast-language:EN-GB">From:</span></b><span lang="EN-US" style="font-size:10.0pt;font-family:"Tahoma","sans-serif";mso-fareast-language:EN-GB">
<a href="mailto:flash-users-bounces@flash.uchicago.edu">flash-users-bounces@flash.uchicago.edu</a> [<a href="mailto:flash-users-bounces@flash.uchicago.edu">mailto:flash-users-bounces@flash.uchicago.edu</a>]
<b>On Behalf Of </b>Jonathan Thurgood<br>
<b>Sent:</b> 04 May 2016 13:16<br>
<b>To:</b> <a href="mailto:flash-users@flash.uchicago.edu">flash-users@flash.uchicago.edu</a><br>
<b>Subject:</b> [FLASH-USERS] USM Custom BC / Grid_bcApplyToRegionSpecialized.F90 / Spurious Current at boundaries for magnetic null point<o:p></o:p></span></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Dear Flash Dev's and Users,<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I've been setting up a USM simulation where the equilibrium is a linear 2D magnetic X-point of the form B=[x,-y,0]) (e.g.,
<a href="http://ukads.nottingham.ac.uk/abs/2011A%26A...533A..18M">http://ukads.nottingham.ac.uk/abs/2011A%26A...533A..18M</a> ,
<a href="http://ukads.nottingham.ac.uk/abs/2004A%26A...420.1129M">http://ukads.nottingham.ac.uk/abs/2004A%26A...420.1129M</a>).
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">The default boundary options don't play well with the null due to the non-constant / curved field through the boundary. The various specifications end up implying some spurious current on the boundaries, and so launch waves into the domain.
To overcome this I need to implement my own boundary conditions in Grid_bcApplyToRegionSpecialised.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">One particular approach I have tried is to force the B-field variables to be their equilibrium values on the boundary (from the analytical prescription), which I believe should properly maintain the equilibrium. Specifically, this boundary
condition is intended to set velocity to zero, other UNK variables similar to "reflecting" (without the sign change) and the B field to B=[x,-y,0].
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">This doesn't seem to be working as I expected - after only coding up the custom BC for the lower face (but not the high faces), I found that that the result is that the current along each boundary is the same at a given output (as opposed
to the lower left corner being different to the upper right, say). This suggests that my "implementation" has some silly mistake and key lines are not being executed when I think they should be.
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I further tested this by going to the lines I had expected to modify the B-field variables and setting the variable to be a large value (1e6) to see if it changes (breaks) the simulation. Many of the lines actually appear to do nothing
/ not be executed and I am not sure why, likely I have not properly understood the logic used in their execution. I have tagged these "!*1e6 does nothing" or "!* 1e6 breaks simulation" below.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">If anyone could give me some pointers with this I would be extremely grateful.
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Best regards,<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Jonathan <o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">(Post-doc working in magnetic reconnection, Northumbria University Solar Physics Group, UK)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">subroutine Grid_bcApplyToRegionSpecialized(bcType,gridDataStruct,&<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> guard,axis,face,regionData,regionSize,mask,applied,&<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> blockHandle,secondDir,thirdDir,endPoints,blkLimitsGC, idest)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">use Grid_data, ONLY : gr_meshMe<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#include "constants.h"<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> implicit none<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> integer, intent(IN) :: bcType,axis,face,guard,gridDataStruct<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> integer,dimension(REGION_DIM),intent(IN) :: regionSize<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> real,dimension(regionSize(BC_DIR),&<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> regionSize(SECOND_DIR),&<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> regionSize(THIRD_DIR),&<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> regionSize(STRUCTSIZE)),intent(INOUT)::regionData<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> logical,intent(IN),dimension(regionSize(STRUCTSIZE)):: mask<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> logical, intent(OUT) :: applied<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> integer,intent(IN) :: blockHandle<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> integer,intent(IN) :: secondDir,thirdDir<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> integer,intent(IN),dimension(LOW:HIGH,MDIM) :: endPoints, blkLimitsGC<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> integer,intent(IN),OPTIONAL:: idest<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">!!!jonathan’s added declarations<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> integer :: i,j, k,ivar,je,ke,n,varCount,bcTypeActual<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> logical :: isFace<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> integer :: sizeX, sizeY, sizeZ<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> real, allocatable,dimension(:) :: xCoord, yCoord, XCoordF, yCoordF<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">! stub provided<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> applied = .false.<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">!added lines to the stub below this comment<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> select case (bcType)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> case(USER_DEFINED)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> applied = .true.<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> case default<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> applied = .false. <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> return <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> end select<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> je=regionSize(SECOND_DIR)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> ke=regionSize(THIRD_DIR)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> varCount=regionSize(STRUCTSIZE)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> isFace = (gridDataStruct==FACEX).and.(axis==IAXIS)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> isFace = isFace.or.((gridDataStruct==FACEY).and.(axis==JAXIS))<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> isFace = isFace.or.((gridDataStruct==FACEZ).and.(axis==KAXIS))<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> !some allocation and calls to populate coordinate arrays to allow analytical prescription of B on boundary<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> sizeX = blkLimitsGC(HIGH,IAXIS)-blkLimitsGC(LOW,IAXIS)+1 !double check sizes<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> sizeY = blkLimitsGC(HIGH,JAXIS)-blkLimitsGC(LOW,JAXIS)+1<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> sizeZ = blkLimitsGC(HIGH,KAXIS)-blkLimitsGC(LOW,KAXIS)+1<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> allocate(xCoord(sizeX) )<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> allocate(yCoord(sizeY) )<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> allocate(xCoordF(sizeX+1) )<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> allocate(yCoordF(sizeY+1) )<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> xCoord = 0.0<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> yCoord = 0.0<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> xCoordF = 0.0<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> yCoordF = 0.0<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> call gr_extendedGetCellCoords(IAXIS, blockHandle, gr_meshMe, CENTER , .true. , xCoord, sizeX)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> call gr_extendedGetCellCoords(JAXIS, blockHandle, gr_meshMe, CENTER , .true. , yCoord, sizeY)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> call gr_extendedGetCellCoords(IAXIS, blockHandle, gr_meshMe, FACES, .true. , xCoordF, sizeX+1)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> call gr_extendedGetCellCoords(JAXIS, blockHandle, gr_meshMe, FACES, .true. , yCoordF, sizeY+1)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> !main loop over variables<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> do ivar = 1,varCount<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if(mask(ivar)) then<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> call gr_bcMapBcType(bcTypeActual,bcType,ivar,gridDataStruct,axis,face,idest)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if(face==LOW) then<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> select case (bcTypeActual)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> case(USER_DEFINED)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> k = 2*guard+1<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if(isFace)k=k+1<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> do i = 1,guard<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> regionData(i,1:je,1:ke,ivar) = regionData(k-i,1:je,1:ke,ivar) ! do znd-type for all variables first, then overwrite
<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> ! e.g. b field with more specific values<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> !*1e6 breaks simulation<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef VELX_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if ( ivar==VELX_VAR ) regionData(i,1:je,1:ke,ivar) = 0.0 !overwrite previously assigned value with zero for velocities !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef VELY_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if ( ivar==VELY_VAR ) regionData(i,1:je,1:ke,ivar) = 0.0 !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef VELZ_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if ( ivar==VELZ_VAR ) regionData(i,1:je,1:ke,ivar) = 0.0 !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> !set B=[x,-y,0] on lower face to make current free boundary
<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef MAGX_VAR <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if( (ivar==MAGX_VAR).AND.(axis=iaxis) ) regionData(i,1:je,1:ke,ivar) = xCoord(i) !bx = x !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if( (ivar==MAGX_VAR).AND.(axis=jaxis) ) then<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> do j=1,je<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> regionData(i,j,1:ke,ivar) = xCoord(j) !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> enddo<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef MAGY_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if( (ivar==MAGY_VAR).AND.(axis=jaxis) ) regionData(i,1:je,1:ke,ivar) = -yCoord(i) !by = -y !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if( (ivar==MAGY_VAR).AND.(axis=iaxis) ) then<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> do j=1,je<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> regionData(i,j,1:ke,ivar) = -yCoord(j) !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> enddo <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef MAG_FACE_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if( (ivar==MAG_FACE_VAR).AND.(axis=iaxis) ) regionData(i,1:je,1:ke,ivar) = xCoordF(i) !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if( (ivar==MAG_FACE_VAR).AND.(axis=jaxis) ) regionData(i,1:je,1:ke,ivar) = -yCoordF(i) !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> !should also do intermediate step mag field on face if doing a custom BC<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> !http://flash.uchicago.edu/pipermail/flash-users/2012-February/001021.html<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef MAGI_FACE_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if( (ivar==MAGI_FACE_VAR).AND.(axis=iaxis) ) regionData(i,1:je,1:ke,ivar) = xCoordF(i) !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if( (ivar==MAGI_FACE_VAR).AND.(axis=jaxis) ) regionData(i,1:je,1:ke,ivar) = -yCoordF(i) !*1e6 does nothing<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> end do<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> end select<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> else !face==HIGH<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> select case (bcTypeActual)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> case(USER_DEFINED)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> k = 2*guard+1<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if(isFace)k=k+1<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> do i = 1,guard<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> regionData(k-i,1:je,1:ke,ivar)= regionData(i,1:je,1:ke,ivar) !znd for all first<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef VELX_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if ( ivar==VELX_VAR ) regionData(k-i,1:je,1:ke,ivar) = 0.0 !over-write with zero for velocities<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef VELY_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if ( ivar==VELY_VAR ) regionData(k-i,1:je,1:ke,ivar) = 0.0 !over-write with zero for velocities<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#ifdef VELZ_VAR<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> if ( ivar==VELZ_VAR ) regionData(k-i,1:je,1:ke,ivar) = 0.0 !over-write with zero for velocities<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">#endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> end do<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> end select<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> endif<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> endif <o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> enddo !end ivar do loop<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> deallocate(xCoord)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> deallocate(yCoord)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> deallocate(xCoordF)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> deallocate(yCoordF)<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt"> return<o:p></o:p></p>
<p class="MsoNormal" style="margin-left:36.0pt">end subroutine Grid_bcApplyToRegionSpecialized<o:p></o:p></p>
</div>
</body>
</html>